From e4da2bbd16bcefc1c85932e47e697e4170870d33 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Tue, 18 Aug 2026 09:55:23 +0100 Subject: [PATCH] Modeling - Modernize BRepGProp_Gauss integration storage (#1473) - Replace handled math_Vector buffers with NCollection_Array1 and NCollection_LocalArray storage. - Reuse Gauss and knot buffers across integration iterations. - Use zero-based size_t indexing and standard algorithms for interval construction and error selection. - Convert the Gauss integration type to enum class. - Simplify inertia initialization and remove obsolete helper methods. --- .../TKTopAlgo/BRepGProp/BRepGProp_Gauss.cxx | 612 ++++++++---------- .../TKTopAlgo/BRepGProp/BRepGProp_Gauss.hxx | 82 +-- .../TKTopAlgo/BRepGProp/BRepGProp_Sinert.cxx | 6 +- .../TKTopAlgo/BRepGProp/BRepGProp_Vinert.cxx | 18 +- .../TKTopAlgo/GTests/BRepGProp_Test.cxx | 236 +++++++ 5 files changed, 559 insertions(+), 395 deletions(-) diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.cxx index 853f746f75..064858d957 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.cxx @@ -11,14 +11,17 @@ // 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 -// If the following is defined the error of algorithm is calculated by static moments -#define IS_MIN_DIM +#include namespace { @@ -28,9 +31,42 @@ static const double EPS_DIM = 1.e-30; static const double ERROR_ALGEBR_RATIO = 2.0 / 3.0; // Maximum of GaussPoints on a subinterval and maximum of subintervals -static const int GPM = math::GaussPointsMax(); -static const int SUBS_POWER = 32; -static const int SM = SUBS_POWER * GPM + 1; +static const int GPM = math::GaussPointsMax(); +static const int SUBS_POWER = 32; +static const size_t SM = static_cast(SUBS_POWER * GPM + 1); + +static constexpr int THE_GAUSS_INLINE_SIZE = 64; +static constexpr int THE_KNOT_INLINE_SIZE = 2; + +template +static NCollection_Array1 prepareArray( + NCollection_LocalArray& theStorage, + const size_t theSize) +{ + if (theStorage.Size() < theSize) + { + theStorage.Reallocate(theSize); + } + return NCollection_Array1(theStorage.begin(), theSize); +} + +static void initGaussPoints(const int theNbPoints, + NCollection_LocalArray& theStorage) +{ + NCollection_Array1 aPointArray = + prepareArray(theStorage, static_cast(theNbPoints)); + math_Vector aPoints(aPointArray.Data(), 1, aPointArray.Length()); + math::GaussPoints(theNbPoints, aPoints); +} + +static void initGaussWeights(const int theNbPoints, + NCollection_LocalArray& theStorage) +{ + NCollection_Array1 aWeightArray = + prepareArray(theStorage, static_cast(theNbPoints)); + math_Vector aWeights(aWeightArray.Data(), 1, aWeightArray.Length()); + math::GaussWeights(theNbPoints, aWeights); +} // Auxiliary inner functions to perform arithmetic operations. static double Add(const double theA, const double theB) @@ -153,148 +189,66 @@ static double MultInf(const double theA, const double theB) return theA * theB; } + +//================================================================================================= + +static size_t maxSubs(const size_t theN, const size_t theCoeff = 32) +{ + return static_cast(IntegerLast()) / theCoeff < theN ? static_cast(IntegerLast()) + : theN * theCoeff + 1; +} + +//================================================================================================= + +template +static size_t fillIntervalBounds(const double theA, + const double theB, + const NCollection_Array1& theKnots, + const size_t theNumSubs, + NCollection_Array1& theInerts, + NCollection_Array1& theParam1, + NCollection_Array1& theParam2, + NCollection_Array1& theError, + NCollection_Array1* theCommonError) +{ + const size_t aSize = std::max(theKnots.Size(), maxSubs(theKnots.Size() - 1, theNumSubs)); + + if (aSize > theParam1.Size()) + { + theInerts.Resize(aSize, false); + theParam1.Resize(aSize, false); + theParam2.Resize(aSize, false); + theError.Resize(aSize, false); + theError.Init(0.0); + + if (theCommonError != nullptr) + { + theCommonError->Resize(aSize, false); + theCommonError->Init(0.0); + } + } + + const auto aFirstKnot = std::upper_bound(theKnots.cbegin(), theKnots.cend(), theA); + const auto aLastKnot = std::lower_bound(aFirstKnot, theKnots.cend(), theB); + const auto aNbKnots = static_cast(aLastKnot - aFirstKnot); + + theParam1.ChangeAt(0) = theA; + std::copy(aFirstKnot, aLastKnot, theParam1.begin() + 1); + std::copy(aFirstKnot, aLastKnot, theParam2.begin()); + theParam2.ChangeAt(aNbKnots) = theB; + return aNbKnots + 1; +} } // namespace -//================================================================================================= - -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) +BRepGProp_Gauss::BRepGProp_Gauss(const GaussType theType) + : myType(theType), + add(::Add), + mult(::Mult) { } //================================================================================================= -void BRepGProp_Gauss::Inertia::Reset() -{ - memset(reinterpret_cast(this), 0, sizeof(BRepGProp_Gauss::Inertia)); -} - -//================================================================================================= - -BRepGProp_Gauss::BRepGProp_Gauss(const BRepGProp_GaussType theType) - : myType(theType) -{ - add = (::Add); - mult = (::Mult); -} - -//================================================================================================= - -int BRepGProp_Gauss::MaxSubs(const int theN, const int theCoeff) -{ - return IntegerLast() / theCoeff < theN ? IntegerLast() : theN * theCoeff + 1; -} - -//================================================================================================= - -void BRepGProp_Gauss::Init(NCollection_Handle& theOutVec, - const double theValue, - const int theFirst, - const int theLast) -{ - if (theLast - theFirst == 0) - { - theOutVec->Init(theValue); - } - else - { - for (int i = theFirst; i <= theLast; ++i) - { - theOutVec->Value(i) = theValue; - } - } -} - -//================================================================================================= - -void BRepGProp_Gauss::InitMass(const double theValue, - const int theFirst, - const int theLast, - InertiaArray& theArray) -{ - if (theArray.IsNull()) - { - return; - } - - int aFirst = theFirst; - int aLast = theLast; - - if (theLast - theFirst == 0) - { - aFirst = theArray->Lower(); - aLast = theArray->Upper(); - } - - for (int i = aFirst; i <= aLast; ++i) - { - theArray->ChangeValue(i).Mass = theValue; - } -} - -//================================================================================================= - -int BRepGProp_Gauss::FillIntervalBounds(const double theA, - const double theB, - const NCollection_Array1& theKnots, - const int theNumSubs, - InertiaArray& theInerts, - NCollection_Handle& theParam1, - NCollection_Handle& theParam2, - NCollection_Handle& theError, - NCollection_Handle& theCommonError) -{ - const int aSize = std::max(theKnots.Upper(), MaxSubs(theKnots.Upper() - 1, theNumSubs)); - - if (aSize - 1 > theParam1->Upper()) - { - theInerts = new NCollection_Array1(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()) - { - theCommonError = new math_Vector(1, aSize, 0.0); - } - } - - int j = 1, k = 1; - theParam1->Value(j++) = theA; - - const int aLength = theKnots.Upper(); - for (int i = 1; i <= aLength; ++i) - { - const double 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; -} - -//================================================================================================= - void BRepGProp_Gauss::computeVInertiaOfElementaryPart(const gp_Pnt& thePoint, const gp_Vec& theNormal, const gp_Pnt& theLocation, @@ -545,36 +499,37 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, double anEpsilon = std::abs(theEps); - BRepGProp_Gauss::Inertia anInertia; - InertiaArray anInertiaL = new NCollection_Array1(1, SM); - InertiaArray anInertiaU = new NCollection_Array1(1, SM); + BRepGProp_Gauss::Inertia anInertia; + NCollection_Array1 anInertiaL(SM); + NCollection_Array1 anInertiaU(SM); // Prepare Gauss points and weights - NCollection_Handle LGaussP[2]; - NCollection_Handle LGaussW[2]; - NCollection_Handle UGaussP[2]; - NCollection_Handle UGaussW[2]; - const int aNbGaussPoint = RealToInt(std::ceil(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); + NCollection_LocalArray LGaussP[2] = { + NCollection_LocalArray(static_cast(GPM)), + NCollection_LocalArray(static_cast(aNbGaussPoint))}; + NCollection_LocalArray LGaussW[2] = { + NCollection_LocalArray(static_cast(GPM)), + NCollection_LocalArray(static_cast(aNbGaussPoint))}; + NCollection_LocalArray UGaussP[2] = { + NCollection_LocalArray(static_cast(GPM)), + NCollection_LocalArray(static_cast(aNbGaussPoint))}; + NCollection_LocalArray UGaussW[2] = { + NCollection_LocalArray(static_cast(GPM)), + NCollection_LocalArray(static_cast(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); + NCollection_Array1 L1(SM); + NCollection_Array1 L2(SM); + NCollection_Array1 U1(SM); + NCollection_Array1 U2(SM); - NCollection_Handle L1 = new math_Vector(1, SM); - NCollection_Handle L2 = new math_Vector(1, SM); - NCollection_Handle U1 = new math_Vector(1, SM); - NCollection_Handle U2 = new math_Vector(1, SM); - - NCollection_Handle ErrL = new math_Vector(1, SM, 0.0); - NCollection_Handle ErrU = new math_Vector(1, SM, 0.0); - NCollection_Handle ErrUL = new math_Vector(1, SM, 0.0); + NCollection_Array1 ErrL(SM); + NCollection_Array1 ErrU(SM); + NCollection_Array1 ErrUL(SM); + ErrL.Init(0.0); + ErrU.Init(0.0); + ErrUL.Init(0.0); // Face parametrization in U and V direction double BV1, BV2, BU1, BU2; @@ -582,7 +537,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, checkBounds(BU1, BU2, BV1, BV2); // - const int NumSubs = SUBS_POWER; + const size_t NumSubs = SUBS_POWER; const TopoDS_Face& aF = theSurface.GetFace(); // clang-format off const bool isNaturalRestriction = (aF.NbChildren () == 0); //theSurface.NaturalRestriction(); @@ -599,9 +554,12 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, gp_Vec2d Vuv; double Dul; // Dul = Du / Dl - int iLS, iLSubEnd, iGL, iGLEnd, NbLGaussP[2], LRange[2], iL, kL, kLEnd, IL, JL; - int i, iUSubEnd, NbUGaussP[2], URange[2], kU, kUEnd, IU, JU; - int UMaxSubs, LMaxSubs; + size_t iLS, iLSubEnd, LRange[2], iL, kL, IL, JL; + size_t iUSubEnd, URange[2], kU, IU, JU; + size_t UMaxSubs, LMaxSubs; + int iGL, iGLEnd, NbLGaussP[2]; + int NbUGaussP[2]; + size_t kLEnd, kUEnd; double ErrorU, ErrorL, ErrorLMax = 0.0, Eps = 0.0, EpsL = 0.0, EpsU = 0.0; iGLEnd = isErrorCalculation ? 2 : 1; @@ -609,15 +567,19 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, NbUGaussP[0] = theSurface.SIntOrder(anEpsilon); NbUGaussP[1] = RealToInt(std::ceil(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]); + initGaussPoints(NbUGaussP[0], UGaussP[0]); + initGaussWeights(NbUGaussP[0], UGaussW[0]); + initGaussPoints(NbUGaussP[1], UGaussP[1]); + initGaussWeights(NbUGaussP[1], UGaussW[1]); - const int aNbUSubs = theSurface.SUIntSubs(); - NCollection_Array1 UKnots(1, aNbUSubs + 1); + const int aNbUSubs = theSurface.SUIntSubs(); + NCollection_LocalArray aUKnotStorage( + static_cast(aNbUSubs + 1)); + NCollection_Array1 UKnots = aUKnotStorage.ToArray1(); + UKnots.UpdateLowerBound(1); theSurface.UKnots(UKnots); + NCollection_LocalArray aLKnotStorage; while (isNaturalRestriction || theDomain.More()) { if (isNaturalRestriction) @@ -635,13 +597,15 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, NbLGaussP[1] = RealToInt(std::ceil(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]); + initGaussPoints(NbLGaussP[0], LGaussP[0]); + initGaussWeights(NbLGaussP[0], LGaussW[0]); + initGaussPoints(NbLGaussP[1], LGaussP[1]); + initGaussWeights(NbLGaussP[1], LGaussW[1]); const int aNbLSubs = isNaturalRestriction ? theSurface.SVIntSubs() : theSurface.LIntSubs(); - NCollection_Array1 LKnots(1, aNbLSubs + 1); + NCollection_Array1 LKnots = + prepareArray(aLKnotStorage, static_cast(aNbLSubs + 1)); + LKnots.UpdateLowerBound(1); if (isNaturalRestriction) { @@ -661,39 +625,42 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, if (std::abs(l2 - l1) > EPS_PARAM) { - iLSubEnd = FillIntervalBounds(l1, l2, LKnots, NumSubs, anInertiaL, L1, L2, ErrL, ErrUL); - LMaxSubs = BRepGProp_Gauss::MaxSubs(iLSubEnd); + iLSubEnd = fillIntervalBounds(l1, l2, LKnots, NumSubs, anInertiaL, L1, L2, ErrL, &ErrUL); + LMaxSubs = 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); + std::for_each(anInertiaL.begin(), anInertiaL.end(), [](Inertia& theInertia) { + theInertia.Mass = 0.0; + }); + ErrL.Init(0.0); + ErrUL.Init(0.0); 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); + LRange[0] = IL = + static_cast(std::max_element(ErrL.cbegin(), ErrL.cend()) - ErrL.cbegin()); + LRange[1] = JL - 1; + L1.ChangeAt(JL - 1) = (L1.At(IL) + L2.At(IL)) * 0.5; + L2.ChangeAt(JL - 1) = L2.At(IL); + L2.ChangeAt(IL) = L1.At(JL - 1); } else { - LRange[0] = IL = JL; + LRange[0] = IL = JL - 1; } - if (JL == LMaxSubs || std::abs(L2->Value(JL) - L1->Value(JL)) < EPS_PARAM) + if (JL == LMaxSubs || std::abs(L2.At(JL - 1) - L1.At(JL - 1)) < EPS_PARAM) { if (kLEnd == 1) { - anInertiaL->ChangeValue(JL).Reset(); - ErrL->Value(JL) = 0.0; + anInertiaL.ChangeAt(JL - 1) = Inertia(); + ErrL.ChangeAt(JL - 1) = 0.0; } else { @@ -708,8 +675,8 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, 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)); + lm = 0.5 * (L2.At(iLS) + L1.At(iLS)); + lr = 0.5 * (L2.At(iLS) - L1.At(iLS)); CIx = CIy = CIz = CIxy = CIxz = CIyz = 0.0; @@ -717,19 +684,19 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, { CDim[iGL] = CIxx[iGL] = CIyy[iGL] = CIzz[iGL] = 0.0; - for (iL = 1; iL <= NbLGaussP[iGL]; iL++) + for (iL = 0; iL < static_cast(NbLGaussP[iGL]); ++iL) { - l = lm + lr * LGaussP[iGL]->Value(iL); + l = lm + lr * LGaussP[iGL].begin()[iL]; if (isNaturalRestriction) { v = l; u2 = BU2; - Dul = LGaussW[iGL]->Value(iL); + Dul = LGaussW[iGL].begin()[iL]; } else { theSurface.D12d(l, Puv, Vuv); - Dul = Vuv.Y() * LGaussW[iGL]->Value(iL); // Dul = Du / Dl + Dul = Vuv.Y() * LGaussW[iGL].begin()[iL]; // Dul = Du / Dl if (std::abs(Dul) < EPS_PARAM) { @@ -759,51 +726,53 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, } } - ErrUL->Value(iLS) = 0.0; - kUEnd = 1; - JU = 0; + ErrUL.ChangeAt(iLS) = 0.0; + kUEnd = 1; + JU = 0; if (std::abs(u2 - u1) < EPS_PARAM) { continue; } - NCollection_Handle aDummy; iUSubEnd = - FillIntervalBounds(u1, u2, UKnots, NumSubs, anInertiaU, U1, U2, ErrU, aDummy); - UMaxSubs = BRepGProp_Gauss::MaxSubs(iUSubEnd); + fillIntervalBounds(u1, u2, UKnots, NumSubs, anInertiaU, U1, U2, ErrU, nullptr); + UMaxSubs = maxSubs(iUSubEnd); if (UMaxSubs > SM) { UMaxSubs = SM; } - BRepGProp_Gauss::InitMass(0.0, 1, UMaxSubs, anInertiaU); - BRepGProp_Gauss::Init(ErrU, 0.0, 1, UMaxSubs); + std::for_each(anInertiaU.begin(), anInertiaU.end(), [](Inertia& theInertia) { + theInertia.Mass = 0.0; + }); + ErrU.Init(0.0); ErrorU = 0.0; do { // while: U if (++JU > iUSubEnd) { - URange[0] = IU = ErrU->Max(); - URange[1] = JU; + URange[0] = IU = static_cast( + std::max_element(ErrU.cbegin(), ErrU.cend()) - ErrU.cbegin()); + URange[1] = JU - 1; - U1->Value(JU) = (U1->Value(IU) + U2->Value(IU)) * 0.5; - U2->Value(JU) = U2->Value(IU); - U2->Value(IU) = U1->Value(JU); + U1.ChangeAt(JU - 1) = (U1.At(IU) + U2.At(IU)) * 0.5; + U2.ChangeAt(JU - 1) = U2.At(IU); + U2.ChangeAt(IU) = U1.At(JU - 1); } else { - URange[0] = IU = JU; + URange[0] = IU = JU - 1; } - if (JU == UMaxSubs || std::abs(U2->Value(JU) - U1->Value(JU)) < EPS_PARAM) + if (JU == UMaxSubs || std::abs(U2.At(JU - 1) - U1.At(JU - 1)) < EPS_PARAM) { if (kUEnd == 1) { - ErrU->Value(JU) = 0.0; - anInertiaU->ChangeValue(JU).Reset(); + ErrU.ChangeAt(JU - 1) = 0.0; + anInertiaU.ChangeAt(JU - 1) = Inertia(); } else { @@ -823,22 +792,22 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, { BRepGProp_Gauss::Inertia aLocal[2]; - int iUS = URange[kU]; - const int aLength = iGLEnd - iGL; + const size_t iUS = URange[kU]; + const int aLength = iGLEnd - iGL; - const double um = 0.5 * (U2->Value(iUS) + U1->Value(iUS)); - const double ur = 0.5 * (U2->Value(iUS) - U1->Value(iUS)); + const double um = 0.5 * (U2.At(iUS) + U1.At(iUS)); + const double ur = 0.5 * (U2.At(iUS) - U1.At(iUS)); - for (int iGU = 0; iGU < aLength; ++iGU) + for (size_t iGU = 0; iGU < static_cast(aLength); ++iGU) { - for (int iU = 1; iU <= NbUGaussP[iGU]; ++iU) + for (size_t iU = 0; iU < static_cast(NbUGaussP[iGU]); ++iU) { - double w = UGaussW[iGU]->Value(iU); - const double u = um + ur * UGaussP[iGU]->Value(iU); + double w = UGaussW[iGU].begin()[iU]; + const double u = um + ur * UGaussP[iGU].begin()[iU]; theSurface.Normal(u, v, aPoint, aNormal); - if (myType == Vinert) + if (myType == GaussType::Vinert) { computeVInertiaOfElementaryPart(aPoint, aNormal, @@ -866,11 +835,11 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, } } - BRepGProp_Gauss::Inertia& anUI = anInertiaU->ChangeValue(iUS); + BRepGProp_Gauss::Inertia& anUI = anInertiaU.ChangeAt(iUS); anUI.Mass = mult(aLocal[0].Mass, ur); - if (myType == Vinert) + if (myType == GaussType::Vinert) { anUI.Ixx = mult(aLocal[0].Ixx, ur); anUI.Iyy = mult(aLocal[0].Iyy, ur); @@ -884,7 +853,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, double aDMass = std::abs(aLocal[1].Mass - aLocal[0].Mass); - if (myType == Vinert) + if (myType == GaussType::Vinert) { aLocal[1].Ixx = std::abs(aLocal[1].Ixx - aLocal[0].Ixx); aLocal[1].Iyy = std::abs(aLocal[1].Iyy - aLocal[0].Iyy); @@ -898,11 +867,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, 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); + ErrU.ChangeAt(iUS) = mult(aDMass, ur); } else { @@ -916,7 +881,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, anUI.Ixz = mult(aLocal[0].Ixz, ur); anUI.Iyz = mult(aLocal[0].Iyz, ur); - ErrU->Value(iUS) = mult(aDMass, ur); + ErrU.ChangeAt(iUS) = mult(aDMass, ur); } } } @@ -924,13 +889,13 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, if (JU == iUSubEnd) { kUEnd = 2; - ErrorU = ErrU->Value(ErrU->Max()); + ErrorU = *std::max_element(ErrU.cbegin(), ErrU.cend()); } } while ((ErrorU - EpsU > 0.0 && EpsU != 0.0) || kUEnd == 1); - for (i = 1; i <= JU; ++i) + for (size_t anIndex = 0; anIndex < JU; ++anIndex) { - const BRepGProp_Gauss::Inertia& anIU = anInertiaU->Value(i); + const BRepGProp_Gauss::Inertia& anIU = anInertiaU.At(anIndex); CDim[iGL] = add(CDim[iGL], mult(anIU.Mass, Dul)); CIxx[iGL] = add(CIxx[iGL], mult(anIU.Ixx, Dul)); @@ -943,11 +908,11 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, continue; } - ErrUL->Value(iLS) = ErrorU * std::abs((u2 - u1) * Dul); + ErrUL.ChangeAt(iLS) = ErrorU * std::abs((u2 - u1) * Dul); - for (i = 1; i <= JU; ++i) + for (size_t anIndex = 0; anIndex < JU; ++anIndex) { - const BRepGProp_Gauss::Inertia& anIU = anInertiaU->Value(i); + const BRepGProp_Gauss::Inertia& anIU = anInertiaU.At(anIndex); CIx = add(CIx, mult(anIU.Ix, Dul)); CIy = add(CIy, mult(anIU.Iy, Dul)); @@ -960,7 +925,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, } // for: iL } // for: iGL - BRepGProp_Gauss::Inertia& aLI = anInertiaL->ChangeValue(iLS); + BRepGProp_Gauss::Inertia& aLI = anInertiaL.ChangeAt(iLS); aLI.Mass = mult(CDim[0], lr); aLI.Ixx = mult(CIxx[0], lr); @@ -971,23 +936,19 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, { double aSubDim = std::abs(CDim[1] - CDim[0]); - if (myType == Vinert) + if (myType == GaussType::Vinert) { - ErrorU = ErrUL->Value(iLS); + ErrorU = ErrUL.At(iLS); CIxx[1] = std::abs(CIxx[1] - CIxx[0]); CIyy[1] = std::abs(CIyy[1] - CIyy[0]); CIzz[1] = std::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); + ErrL.ChangeAt(iLS) = add(mult(aSubDim, lr), ErrorU); } else { - ErrL->Value(iLS) = add(mult(aSubDim, lr), ErrUL->Value(iLS)); + ErrL.ChangeAt(iLS) = add(mult(aSubDim, lr), ErrUL.At(iLS)); } } @@ -1008,30 +969,11 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, kLEnd = 2; double DDim = 0.0; - for (i = 1; i <= JL; ++i) + for (size_t anIndex = 0; anIndex < JL; ++anIndex) { - DDim += anInertiaL->Value(i).Mass; + DDim += anInertiaL.At(anIndex).Mass; } -#ifndef IS_MIN_DIM - { - if (myType == Vinert) - { - double 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 = std::abs(DIxx) + std::abs(DIyy) + std::abs(DIzz); - } - } -#endif - DDim = std::abs(DDim * anEpsilon); if (DDim > Eps) @@ -1042,13 +984,13 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, } if (kLEnd == 2) { - ErrorL = ErrL->Value(ErrL->Max()); + ErrorL = *std::max_element(ErrL.cbegin(), ErrL.cend()); } } while ((ErrorL - EpsL > 0.0 && isVerifyComputation) || kLEnd == 1); - for (i = 1; i <= JL; i++) + for (size_t anIndex = 0; anIndex < JL; ++anIndex) { - addAndRestoreInertia(anInertiaL->Value(i), anInertia); + addAndRestoreInertia(anInertiaL.At(anIndex), anInertia); } ErrorLMax = std::max(ErrorLMax, ErrorL); @@ -1062,7 +1004,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, theDomain.Next(); } - if (myType == Vinert) + if (myType == GaussType::Vinert) { convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass); } @@ -1076,14 +1018,6 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, if (theOutMass != 0.0) { Eps = ErrorLMax / std::abs(theOutMass); - -#ifndef IS_MIN_DIM - { - if (myType == Vinert) - Eps = ErrorLMax - / (std::abs(anInertia.Ixx) + std::abs(anInertia.Iyy) + std::abs(anInertia.Izz)); - } -#endif } else { @@ -1108,7 +1042,7 @@ double BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, gp_Pnt& theOutGravityCenter, gp_Mat& theOutInertia) { - Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type"); + Standard_ASSERT_RAISE(myType == GaussType::Sinert, "BRepGProp_Gauss: Incorrect type"); return Compute(theSurface, theDomain, @@ -1130,7 +1064,7 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, gp_Pnt& theOutGravityCenter, gp_Mat& theOutInertia) { - Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type"); + Standard_ASSERT_RAISE(myType == GaussType::Sinert, "BRepGProp_Gauss: Incorrect type"); double u1, u2, v1, v2; theSurface.Bounds(u1, u2, v1, v2); @@ -1143,12 +1077,18 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, const int NbGaussgp_Pnts = std::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); + NCollection_LocalArray aGaussSPVStorage( + static_cast(NbGaussgp_Pnts)); + NCollection_LocalArray aGaussSWVStorage( + static_cast(NbGaussgp_Pnts)); + NCollection_Array1 GaussSPV = aGaussSPVStorage.ToArray1(); + NCollection_Array1 GaussSWV = aGaussSWVStorage.ToArray1(); + initGaussPoints(NbGaussgp_Pnts, aGaussSPVStorage); + initGaussWeights(NbGaussgp_Pnts, aGaussSWVStorage); - BRepGProp_Gauss::Inertia anInertia; + BRepGProp_Gauss::Inertia anInertia; + NCollection_LocalArray aGaussCPStorage; + NCollection_LocalArray aGaussCWStorage; for (; theDomain.More(); theDomain.Next()) { if (!theSurface.Load(theDomain.Value())) @@ -1160,10 +1100,12 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, NbCGaussgp_Pnts = std::max(NbCGaussgp_Pnts, NbGaussgp_Pnts); - math_Vector GaussCP(1, NbCGaussgp_Pnts); - math_Vector GaussCW(1, NbCGaussgp_Pnts); - math::GaussPoints(NbCGaussgp_Pnts, GaussCP); - math::GaussWeights(NbCGaussgp_Pnts, GaussCW); + initGaussPoints(NbCGaussgp_Pnts, aGaussCPStorage); + initGaussWeights(NbCGaussgp_Pnts, aGaussCWStorage); + NCollection_Array1 GaussCP = + prepareArray(aGaussCPStorage, static_cast(NbCGaussgp_Pnts)); + NCollection_Array1 GaussCW = + prepareArray(aGaussCWStorage, static_cast(NbCGaussgp_Pnts)); const double l1 = theSurface.FirstParameter(); const double l2 = theSurface.LastParameter(); @@ -1171,9 +1113,9 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, const double lr = 0.5 * (l2 - l1); BRepGProp_Gauss::Inertia aCInertia; - for (int i = 1; i <= NbCGaussgp_Pnts; ++i) + for (size_t i = 0; i < static_cast(NbCGaussgp_Pnts); ++i) { - const double l = lm + lr * GaussCP(i); + const double l = lm + lr * GaussCP.At(i); gp_Pnt2d Puv; gp_Vec2d Vuv; @@ -1182,15 +1124,15 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, const double v = Puv.Y(); u2 = Puv.X(); - const double Dul = Vuv.Y() * GaussCW(i); + const double Dul = Vuv.Y() * GaussCW.At(i); const double um = 0.5 * (u2 + u1); const double ur = 0.5 * (u2 - u1); BRepGProp_Gauss::Inertia aLocalInertia; - for (int j = 1; j <= NbGaussgp_Pnts; ++j) + for (size_t j = 0; j < static_cast(NbGaussgp_Pnts); ++j) { - const double u = add(um, mult(ur, GaussSPV(j))); - const double aWeight = Dul * GaussSWV(j); + const double u = add(um, mult(ur, GaussSPV.At(j))); + const double aWeight = Dul * GaussSWV.At(j); gp_Pnt aPoint; gp_Vec aNormal; @@ -1221,7 +1163,7 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, gp_Pnt& theOutGravityCenter, gp_Mat& theOutInertia) { - Standard_ASSERT_RAISE(myType == Vinert, "BRepGProp_Gauss: Incorrect type"); + Standard_ASSERT_RAISE(myType == GaussType::Vinert, "BRepGProp_Gauss: Incorrect type"); double u1, v1, u2, v2; theSurface.Bounds(u1, u2, v1, v2); @@ -1229,7 +1171,9 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, double _u2 = u2; // OCC104 - BRepGProp_Gauss::Inertia anInertia; + BRepGProp_Gauss::Inertia anInertia; + NCollection_LocalArray aGaussPStorage; + NCollection_LocalArray aGaussWStorage; for (; theDomain.More(); theDomain.Next()) { if (!theSurface.Load(theDomain.Value())) @@ -1242,10 +1186,12 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, const int aNbGaussgp_Pnts = std::min(std::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); + initGaussPoints(aNbGaussgp_Pnts, aGaussPStorage); + initGaussWeights(aNbGaussgp_Pnts, aGaussWStorage); + NCollection_Array1 GaussP = + prepareArray(aGaussPStorage, static_cast(aNbGaussgp_Pnts)); + NCollection_Array1 GaussW = + prepareArray(aGaussWStorage, static_cast(aNbGaussgp_Pnts)); const double l1 = theSurface.FirstParameter(); const double l2 = theSurface.LastParameter(); @@ -1253,9 +1199,9 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, const double lr = 0.5 * (l2 - l1); BRepGProp_Gauss::Inertia aCInertia; - for (int i = 1; i <= aNbGaussgp_Pnts; ++i) + for (size_t i = 0; i < static_cast(aNbGaussgp_Pnts); ++i) { - const double l = lm + lr * GaussP(i); + const double l = lm + lr * GaussP.At(i); gp_Pnt2d Puv; gp_Vec2d Vuv; @@ -1266,15 +1212,15 @@ void BRepGProp_Gauss::Compute(BRepGProp_Face& theSurface, u2 = std::min(std::max(u1, u2), _u2); // OCC104 const double v = std::min(std::max(Puv.Y(), v1), v2); - const double Dul = Vuv.Y() * GaussW(i); + const double Dul = Vuv.Y() * GaussW.At(i); const double um = 0.5 * (u2 + u1); const double ur = 0.5 * (u2 - u1); BRepGProp_Gauss::Inertia aLocalInertia; - for (int j = 1; j <= aNbGaussgp_Pnts; ++j) + for (size_t j = 0; j < static_cast(aNbGaussgp_Pnts); ++j) { - const double u = um + ur * GaussP(j); - const double aWeight = Dul * GaussW(j); + const double u = um + ur * GaussP.At(j); + const double aWeight = Dul * GaussW.At(j); gp_Pnt aPoint; gp_Vec aNormal; @@ -1319,15 +1265,23 @@ void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface, const int VOrder = std::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); + NCollection_LocalArray aGaussPUStorage( + static_cast(UOrder)); + NCollection_LocalArray aGaussWUStorage( + static_cast(UOrder)); + NCollection_LocalArray aGaussPVStorage( + static_cast(VOrder)); + NCollection_LocalArray aGaussWVStorage( + static_cast(VOrder)); + NCollection_Array1 GaussPU = aGaussPUStorage.ToArray1(); + NCollection_Array1 GaussWU = aGaussWUStorage.ToArray1(); + NCollection_Array1 GaussPV = aGaussPVStorage.ToArray1(); + NCollection_Array1 GaussWV = aGaussWVStorage.ToArray1(); - math::GaussPoints(UOrder, GaussPU); - math::GaussWeights(UOrder, GaussWU); - math::GaussPoints(VOrder, GaussPV); - math::GaussWeights(VOrder, GaussWV); + initGaussPoints(UOrder, aGaussPUStorage); + initGaussWeights(UOrder, aGaussWUStorage); + initGaussPoints(VOrder, aGaussPVStorage); + initGaussWeights(VOrder, aGaussWVStorage); const double um = 0.5 * add(UpperU, LowerU); const double vm = 0.5 * add(UpperV, LowerV); @@ -1338,18 +1292,18 @@ void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface, gp_Vec aNormal; BRepGProp_Gauss::Inertia anInertia; - for (int j = 1; j <= VOrder; ++j) + for (size_t j = 0; j < static_cast(VOrder); ++j) { BRepGProp_Gauss::Inertia anInertiaOfElementaryPart; - const double v = add(vm, mult(vr, GaussPV(j))); + const double v = add(vm, mult(vr, GaussPV.At(j))); - for (int i = 1; i <= UOrder; ++i) + for (size_t i = 0; i < static_cast(UOrder); ++i) { - const double aWeight = GaussWU(i); - const double u = add(um, mult(ur, GaussPU(i))); + const double aWeight = GaussWU.At(i); + const double u = add(um, mult(ur, GaussPU.At(i))); theSurface.Normal(u, v, aPoint, aNormal); - if (myType == Vinert) + if (myType == GaussType::Vinert) { computeVInertiaOfElementaryPart(aPoint, aNormal, @@ -1369,7 +1323,7 @@ void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface, } } - multAndRestoreInertia(GaussWV(j), anInertiaOfElementaryPart); + multAndRestoreInertia(GaussWV.At(j), anInertiaOfElementaryPart); addAndRestoreInertia(anInertiaOfElementaryPart, anInertia); } vr = mult(vr, ur); @@ -1380,7 +1334,7 @@ void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface, anInertia.Ixz = mult(vr, anInertia.Ixz); anInertia.Iyz = mult(vr, anInertia.Iyz); - if (myType == Vinert) + if (myType == GaussType::Vinert) { convert(anInertia, theCoeff, theIsByPoint, theOutGravityCenter, theOutInertia, theOutMass); } @@ -1400,7 +1354,7 @@ void BRepGProp_Gauss::Compute(const BRepGProp_Face& theSurface, gp_Pnt& theOutGravityCenter, gp_Mat& theOutInertia) { - Standard_ASSERT_RAISE(myType == Sinert, "BRepGProp_Gauss: Incorrect type"); + Standard_ASSERT_RAISE(myType == GaussType::Sinert, "BRepGProp_Gauss: Incorrect type"); Compute(theSurface, theLocation, nullptr, true, theOutMass, theOutGravityCenter, theOutInertia); } diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.hxx b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.hxx index 75bac0bcb5..43573750b4 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.hxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Gauss.hxx @@ -14,48 +14,44 @@ #ifndef _BRepGProp_Gauss_HeaderFile #define _BRepGProp_Gauss_HeaderFile -#include -#include +#include -template -class math_VectorBase; -using math_Vector = math_VectorBase; +#include + +class BRepGProp_Domain; +class BRepGProp_Face; +class gp_Mat; +class gp_Pnt; +class gp_Vec; //! 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 { + using GaussFunc = std::function; + //! Auxiliary structure for storing of inertial moments. struct Inertia { //! Mass of the current system (without density). //! May correspond to: length, area, volume. - double Mass; + double Mass = 0.0; //! Static moments of inertia. - double Ix; - double Iy; - double Iz; + double Ix = 0.0; + double Iy = 0.0; + double Iz = 0.0; //! Quadratic moments of inertia. - double Ixx; - double Iyy; - double Izz; - double Ixy; - double Ixz; - double Iyz; - - //! Default constructor. - Inertia(); - - //! Zeroes all values. - void Reset(); + double Ixx = 0.0; + double Iyy = 0.0; + double Izz = 0.0; + double Ixy = 0.0; + double Ixz = 0.0; + double Iyz = 0.0; }; - typedef NCollection_Handle> InertiaArray; - typedef double (*BRepGProp_GaussFunc)(const double, const double); - public: //! @name public API //! Describes types of geometric objects. //! - Vinert is 3D closed region of space delimited with: @@ -63,14 +59,14 @@ public: //! @name public API //! -- Point and Surface; //! -- Plane and Surface. //! - Sinert is face in 3D space. - typedef enum + enum class GaussType { Vinert = 0, Sinert - } BRepGProp_GaussType; + }; //! Constructor - Standard_EXPORT explicit BRepGProp_Gauss(const BRepGProp_GaussType theType); + Standard_EXPORT explicit BRepGProp_Gauss(const 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. @@ -233,32 +229,10 @@ private: //! @name private methods gp_Mat& theOutMatrixOfInertia, double& theOutMass); - static int MaxSubs(const int theN, const int theCoeff = 32); - - static void Init(NCollection_Handle& theOutVec, - const double theValue, - const int theFirst = 0, - const int theLast = 0); - - static void InitMass(const double theValue, - const int theFirst, - const int theLast, - InertiaArray& theArray); - - static int FillIntervalBounds(const double theA, - const double theB, - const NCollection_Array1& theKnots, - const int theNumSubs, - InertiaArray& theInerts, - NCollection_Handle& theParam1, - NCollection_Handle& theParam2, - NCollection_Handle& theError, - NCollection_Handle& 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 +private: //! @name private fields + GaussType myType; //!< Type of geometric object + GaussFunc add; //!< Addition operation + GaussFunc mult; //!< Multiplication operation }; -#endif \ No newline at end of file +#endif diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Sinert.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Sinert.cxx index 142d249739..b1c01c3183 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Sinert.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Sinert.cxx @@ -77,7 +77,7 @@ void BRepGProp_Sinert::Perform(const BRepGProp_Face& theSurface) { myEpsilon = 1.0; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Sinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Sinert); aGauss.Compute(theSurface, loc, dim, g, inertia); } @@ -87,7 +87,7 @@ void BRepGProp_Sinert::Perform(BRepGProp_Face& theSurface, BRepGProp_Domain& the { myEpsilon = 1.0; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Sinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Sinert); aGauss.Compute(theSurface, theDomain, loc, dim, g, inertia); } @@ -105,7 +105,7 @@ double BRepGProp_Sinert::Perform(BRepGProp_Face& theSurface, BRepGProp_Domain& theDomain, const double theEps) { - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Sinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Sinert); return myEpsilon = aGauss.Compute(theSurface, theDomain, loc, theEps, dim, g, inertia); } diff --git a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Vinert.cxx b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Vinert.cxx index e8788fe860..116b44f861 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Vinert.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/BRepGProp/BRepGProp_Vinert.cxx @@ -175,7 +175,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, const double theEps) { const double aCoeff[] = {0.0, 0.0, 0.0}; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); return myEpsilon = aGauss.Compute(theSurface, theDomain, loc, theEps, aCoeff, true, dim, g, inertia); @@ -186,7 +186,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface) { const double aCoeff[] = {0.0, 0.0, 0.0}; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); myEpsilon = 1.0; aGauss.Compute(theSurface, loc, aCoeff, true, dim, g, inertia); @@ -197,7 +197,7 @@ void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface) void BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, BRepGProp_Domain& theDomain) { const double aCoeff[] = {0.0, 0.0, 0.0}; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); myEpsilon = 1.0; aGauss.Compute(theSurface, theDomain, loc, aCoeff, true, dim, g, inertia); @@ -224,7 +224,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, theOrigin.Y() - loc.Y(), theOrigin.Z() - loc.Z()}; - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); return myEpsilon = aGauss.Compute(theSurface, theDomain, loc, theEps, aCoeff, true, dim, g, inertia); @@ -234,7 +234,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface, const gp_Pnt& theOrigin) { - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); const double aCoeff[] = {theOrigin.X() - loc.X(), theOrigin.Y() - loc.Y(), theOrigin.Z() - loc.Z()}; @@ -249,7 +249,7 @@ void BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, BRepGProp_Domain& theDomain, const gp_Pnt& theOrigin) { - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); const double aCoeff[] = {theOrigin.X() - loc.X(), theOrigin.Y() - loc.Y(), theOrigin.Z() - loc.Z()}; @@ -279,7 +279,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, 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(); - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); return myEpsilon = aGauss.Compute(theSurface, theDomain, loc, theEps, aCoeff, false, dim, g, inertia); @@ -289,7 +289,7 @@ double BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, void BRepGProp_Vinert::Perform(const BRepGProp_Face& theSurface, const gp_Pln& thePlane) { - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); double aCoeff[4]; thePlane.Coefficients(aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]); @@ -306,7 +306,7 @@ void BRepGProp_Vinert::Perform(BRepGProp_Face& theSurface, BRepGProp_Domain& theDomain, const gp_Pln& thePlane) { - BRepGProp_Gauss aGauss(BRepGProp_Gauss::Vinert); + BRepGProp_Gauss aGauss(BRepGProp_Gauss::GaussType::Vinert); double aCoeff[4]; thePlane.Coefficients(aCoeff[0], aCoeff[1], aCoeff[2], aCoeff[3]); diff --git a/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepGProp_Test.cxx b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepGProp_Test.cxx index 189b286dd9..5b7eb29fa5 100644 --- a/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepGProp_Test.cxx +++ b/src/ModelingAlgorithms/TKTopAlgo/GTests/BRepGProp_Test.cxx @@ -13,21 +13,36 @@ #include #include +#include +#include #include +#include +#include +#include +#include #include #include #include #include +#include #include #include #include +#include #include #include +#include +#include +#include +#include +#include #include #include #include #include +#include #include +#include #include #include #include @@ -41,6 +56,31 @@ #include +namespace +{ +void expectGPropsNear(const GProp_GProps& theActual, + const GProp_GProps& theExpected, + const double theTolerance) +{ + EXPECT_NEAR(theActual.Mass(), theExpected.Mass(), theTolerance); + EXPECT_NEAR(theActual.CentreOfMass().X(), theExpected.CentreOfMass().X(), theTolerance); + EXPECT_NEAR(theActual.CentreOfMass().Y(), theExpected.CentreOfMass().Y(), theTolerance); + EXPECT_NEAR(theActual.CentreOfMass().Z(), theExpected.CentreOfMass().Z(), theTolerance); + + const gp_Mat anActualInertia = theActual.MatrixOfInertia(); + const gp_Mat anExpectedInertia = theExpected.MatrixOfInertia(); + for (int aRow = 1; aRow <= 3; ++aRow) + { + for (int aColumn = 1; aColumn <= 3; ++aColumn) + { + EXPECT_NEAR(anActualInertia.Value(aRow, aColumn), + anExpectedInertia.Value(aRow, aColumn), + theTolerance); + } + } +} +} // namespace + TEST(BRepGPropTest, LinearProperties_EdgeLength) { gp_Pnt aP1(0.0, 0.0, 0.0); @@ -107,6 +147,202 @@ TEST(BRepGPropTest, VolumeProperties_BoxCenterOfMass) EXPECT_NEAR(aCOM.Z(), 5.0, Precision::Confusion()); } +TEST(BRepGPropTest, AdaptiveSurfaceProperties_Box) +{ + const TopoDS_Shape aBox = BRepPrimAPI_MakeBox(10.0, 20.0, 30.0).Shape(); + + GProp_GProps aProps; + BRepGProp::SurfaceProperties(aBox, aProps, 1.0e-6); + + EXPECT_NEAR(aProps.Mass(), 2200.0, 1.0e-9); +} + +TEST(BRepGPropTest, AdaptiveVolumeProperties_Box) +{ + const TopoDS_Shape aBox = BRepPrimAPI_MakeBox(10.0, 20.0, 30.0).Shape(); + + GProp_GProps aProps; + BRepGProp::VolumeProperties(aBox, aProps, 1.0e-6); + + EXPECT_NEAR(aProps.Mass(), 6000.0, 1.0e-9); +} + +TEST(BRepGPropTest, AdaptiveNaturalBSplineSurface_ResizesIntervalArrays) +{ + const int aNbUIntervals = math::GaussPointsMax() + 1; + const int aNbUPoles = aNbUIntervals + 1; + NCollection_Array2 aPoles(1, aNbUPoles, 1, 2); + for (int aUIndex = 1; aUIndex <= aNbUPoles; ++aUIndex) + { + aPoles(aUIndex, 1) = gp_Pnt(static_cast(aUIndex - 1), 0.0, 0.0); + aPoles(aUIndex, 2) = gp_Pnt(static_cast(aUIndex - 1), 2.0, 0.0); + } + + NCollection_Array1 aUKnots(1, aNbUPoles); + NCollection_Array1 aUMultiplicities(1, aNbUPoles); + for (int anIndex = 1; anIndex <= aNbUPoles; ++anIndex) + { + aUKnots(anIndex) = static_cast(anIndex - 1); + aUMultiplicities(anIndex) = anIndex == 1 || anIndex == aNbUPoles ? 2 : 1; + } + + NCollection_Array1 aVKnots(1, 2); + aVKnots(1) = 0.0; + aVKnots(2) = 1.0; + NCollection_Array1 aVMultiplicities(1, 2); + aVMultiplicities.Init(2); + + occ::handle aSurface = + new Geom_BSplineSurface(aPoles, aUKnots, aVKnots, aUMultiplicities, aVMultiplicities, 1, 1); + TopoDS_Face aFace; + BRep_Builder().MakeFace(aFace, aSurface, Precision::Confusion()); + ASSERT_EQ(aFace.NbChildren(), 0); + + BRepGProp_Face aFixedFace(aFace); + EXPECT_EQ(aFixedFace.SUIntSubs(), aNbUIntervals); + BRepGProp_Sinert aFixed(aFixedFace, gp_Pnt(0.0, 0.0, 0.0)); + BRepGProp_Face anAdaptiveFace(aFace); + BRepGProp_Sinert anAdaptive(anAdaptiveFace, gp_Pnt(0.0, 0.0, 0.0), 1.0e-8); + + EXPECT_NEAR(anAdaptive.Mass(), 2.0 * aNbUIntervals, 1.0e-10); + EXPECT_NEAR(anAdaptive.CentreOfMass().X(), 0.5 * aNbUIntervals, 1.0e-10); + EXPECT_NEAR(anAdaptive.CentreOfMass().Y(), 1.0, 1.0e-12); + EXPECT_NEAR(anAdaptive.CentreOfMass().Z(), 0.0, 1.0e-12); + expectGPropsNear(anAdaptive, aFixed, 1.0e-9); +} + +TEST(BRepGPropTest, AdaptiveTrimmedQuarterDisk_ReusesBoundaryBuffers) +{ + const double aRadius = 4.0; + const gp_Pnt anOrigin(0.0, 0.0, 0.0); + const gp_Pnt anArcStart(aRadius, 0.0, 0.0); + const gp_Pnt anArcEnd(0.0, aRadius, 0.0); + + const occ::handle anArc = + GC_MakeArcOfCircle(gp_Circ(gp_Ax2(anOrigin, gp_Dir(0.0, 0.0, 1.0)), aRadius), 0.0, M_PI_2, true) + .Value(); + BRepBuilderAPI_MakeWire aWireBuilder; + aWireBuilder.Add(BRepBuilderAPI_MakeEdge(anArc).Edge()); + aWireBuilder.Add(BRepBuilderAPI_MakeEdge(anArcEnd, anOrigin).Edge()); + aWireBuilder.Add(BRepBuilderAPI_MakeEdge(anOrigin, anArcStart).Edge()); + ASSERT_TRUE(aWireBuilder.IsDone()); + + BRepBuilderAPI_MakeFace aFaceBuilder(gp_Pln(anOrigin, gp_Dir(0.0, 0.0, 1.0)), + aWireBuilder.Wire()); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + + BRepGProp_Face anOrderFace(aFace); + BRepGProp_Domain anOrderDomain(aFace); + int aFirstOrder = 0; + bool hasDifferentOrders = false; + size_t aNbBoundaries = 0; + for (; anOrderDomain.More(); anOrderDomain.Next(), ++aNbBoundaries) + { + ASSERT_TRUE(anOrderFace.Load(anOrderDomain.Value())); + const int anOrder = anOrderFace.IntegrationOrder(); + if (aNbBoundaries == 0) + { + aFirstOrder = anOrder; + } + else + { + hasDifferentOrders = hasDifferentOrders || anOrder != aFirstOrder; + } + } + EXPECT_EQ(aNbBoundaries, 3u); + EXPECT_TRUE(hasDifferentOrders); + + BRepGProp_Face aFixedFace(aFace); + BRepGProp_Domain aFixedDomain(aFace); + BRepGProp_Sinert aFixed(aFixedFace, aFixedDomain, anOrigin); + + BRepGProp_Face anAdaptiveFace(aFace); + BRepGProp_Domain anAdaptiveDomain(aFace); + BRepGProp_Sinert anAdaptive(anAdaptiveFace, anAdaptiveDomain, anOrigin, 1.0e-8); + + EXPECT_GE(anAdaptive.GetEpsilon(), 0.0); + EXPECT_LE(anAdaptive.GetEpsilon(), 1.0e-8); + EXPECT_NEAR(anAdaptive.Mass(), M_PI * aRadius * aRadius / 4.0, 1.0e-10); + EXPECT_NEAR(anAdaptive.CentreOfMass().X(), 4.0 * aRadius / (3.0 * M_PI), 1.0e-7); + EXPECT_NEAR(anAdaptive.CentreOfMass().Y(), 4.0 * aRadius / (3.0 * M_PI), 1.0e-7); + EXPECT_NEAR(anAdaptive.CentreOfMass().Z(), 0.0, 1.0e-12); + expectGPropsNear(anAdaptive, aFixed, 5.0e-5); + + BRepGProp_Face aFallbackFace(aFace); + BRepGProp_Domain aFallbackDomain(aFace); + BRepGProp_Sinert aFallback(aFallbackFace, aFallbackDomain, anOrigin, 1.0e-2); + EXPECT_DOUBLE_EQ(aFallback.GetEpsilon(), 1.0e-2); + expectGPropsNear(aFallback, aFixed, 5.0e-5); +} + +TEST(BRepGPropTest, AdaptiveVolumeProperties_Cylinder) +{ + const double aRadius = 3.0; + const double aHeight = 7.0; + const TopoDS_Shape aCylinder = BRepPrimAPI_MakeCylinder(aRadius, aHeight).Shape(); + + GProp_GProps aFixed; + BRepGProp::VolumeProperties(aCylinder, aFixed); + GProp_GProps anAdaptive; + const double anError = BRepGProp::VolumeProperties(aCylinder, anAdaptive, 1.0e-8); + + EXPECT_GE(anError, 0.0); + EXPECT_LE(anError, 1.0e-8); + EXPECT_NEAR(anAdaptive.Mass(), M_PI * aRadius * aRadius * aHeight, 1.0e-10); + EXPECT_NEAR(anAdaptive.CentreOfMass().X(), 0.0, 1.0e-12); + EXPECT_NEAR(anAdaptive.CentreOfMass().Y(), 0.0, 1.0e-12); + EXPECT_NEAR(anAdaptive.CentreOfMass().Z(), aHeight * 0.5, 1.0e-12); + const double aMass = M_PI * aRadius * aRadius * aHeight; + EXPECT_NEAR(anAdaptive.MatrixOfInertia().Value(1, 1), + aMass * (3.0 * aRadius * aRadius + aHeight * aHeight) / 12.0, + 2.0e-4); + EXPECT_NEAR(anAdaptive.MatrixOfInertia().Value(2, 2), + aMass * (3.0 * aRadius * aRadius + aHeight * aHeight) / 12.0, + 2.0e-4); + EXPECT_NEAR(anAdaptive.MatrixOfInertia().Value(3, 3), aMass * aRadius * aRadius / 2.0, 2.0e-4); + expectGPropsNear(anAdaptive, aFixed, 2.0e-4); +} + +TEST(BRepGPropTest, AdaptivePlaneLimitedInertia) +{ + const gp_Dir aNormal(0.0, 0.0, 1.0); + BRepBuilderAPI_MakeFace aFaceBuilder(gp_Pln(gp_Pnt(0.0, 0.0, 2.0), aNormal), 0.0, 3.0, 0.0, 4.0); + ASSERT_TRUE(aFaceBuilder.IsDone()); + const TopoDS_Face aFace = aFaceBuilder.Face(); + const gp_Pln aReferencePlane(gp_Pnt(0.0, 0.0, 0.0), aNormal); + const gp_Pnt aLocation(0.0, 0.0, 0.0); + + BRepGProp_Face aFixedFace(aFace); + BRepGProp_Domain aFixedDomain(aFace); + BRepGProp_Vinert aFixed(aFixedFace, aFixedDomain, aReferencePlane, aLocation); + + BRepGProp_Face anAdaptiveFace(aFace); + BRepGProp_Domain anAdaptiveDomain(aFace); + BRepGProp_Vinert anAdaptive(anAdaptiveFace, anAdaptiveDomain, aReferencePlane, aLocation, 1.0e-8); + + const gp_Mat anAdaptiveInertia = anAdaptive.MatrixOfInertia(); + EXPECT_NEAR(anAdaptiveInertia.Value(1, 1), 160.0, 1.0e-10); + EXPECT_NEAR(anAdaptiveInertia.Value(2, 2), 104.0, 1.0e-10); + EXPECT_NEAR(anAdaptiveInertia.Value(3, 3), 200.0, 1.0e-10); + EXPECT_NEAR(anAdaptiveInertia.Value(1, 2), -72.0, 1.0e-10); + EXPECT_NEAR(anAdaptiveInertia.Value(1, 3), -36.0, 1.0e-10); + EXPECT_NEAR(anAdaptiveInertia.Value(2, 3), -48.0, 1.0e-10); + expectGPropsNear(anAdaptive, aFixed, 1.0e-10); +} + +TEST(BRepGPropTest, NaturalInfinitePlane_UsesInfiniteArithmetic) +{ + const occ::handle aPlane = + new Geom_Plane(gp_Pnt(0.0, 0.0, 0.0), gp_Dir(0.0, 0.0, 1.0)); + TopoDS_Face aFace; + BRep_Builder().MakeFace(aFace, aPlane, Precision::Confusion()); + + BRepGProp_Face aFaceTool(aFace); + BRepGProp_Sinert aProperties(aFaceTool, gp_Pnt(0.0, 0.0, 0.0)); + EXPECT_TRUE(Precision::IsPositiveInfinite(aProperties.Mass())); +} + TEST(BRepGPropTest, LinearProperties_SkipShared) { BRepPrimAPI_MakeBox aBox(10.0, 10.0, 10.0);