diff --git a/src/ModelingData/TKGeomBase/GTests/FILES.cmake b/src/ModelingData/TKGeomBase/GTests/FILES.cmake index 864bd9fed1..4984b0ba5a 100644 --- a/src/ModelingData/TKGeomBase/GTests/FILES.cmake +++ b/src/ModelingData/TKGeomBase/GTests/FILES.cmake @@ -40,6 +40,7 @@ set(OCCT_TKGeomBase_GTests_FILES GeomConvert_CompCurveToBSplineCurve_Test.cxx Geom2dConvert_CompCurveToBSplineCurve_Test.cxx GeomConvert_Test.cxx + GeomLib_CheckCurveOnSurface_Test.cxx GProp_PEquation_Test.cxx GProp_PGProps_Test.cxx Hermit_Test.cxx diff --git a/src/ModelingData/TKGeomBase/GTests/GeomLib_CheckCurveOnSurface_Test.cxx b/src/ModelingData/TKGeomBase/GTests/GeomLib_CheckCurveOnSurface_Test.cxx new file mode 100644 index 0000000000..d64f872391 --- /dev/null +++ b/src/ModelingData/TKGeomBase/GTests/GeomLib_CheckCurveOnSurface_Test.cxx @@ -0,0 +1,110 @@ +// 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 +#include +#include + +namespace +{ +// Circle of radius theRadius at height theHeight on a cylinder of radius +// theCylRadius around OZ; the pcurve is the corresponding iso line in UV. +GeomLib_CheckCurveOnSurface makeCircleOnCylinderCheck(const double theCircRadius, + const double theCylRadius, + const double theHeight, + double& theMaxDistance) +{ + const double aFirst = 0.0, aLast = 2.0 * M_PI; + + const occ::handle aCirc = new Geom_Circle( + gp_Ax2(gp_Pnt(0.0, 0.0, theHeight), gp_Dir(0.0, 0.0, 1.0), gp_Dir(1.0, 0.0, 0.0)), + theCircRadius); + const occ::handle aC3d = new GeomAdaptor_Curve(aCirc, aFirst, aLast); + + const occ::handle aCyl = new Geom_CylindricalSurface( + gp_Ax3(gp_Pnt(0.0, 0.0, 0.0), gp_Dir(0.0, 0.0, 1.0), gp_Dir(1.0, 0.0, 0.0)), + theCylRadius); + const occ::handle aPLine = + new Geom2d_Line(gp_Pnt2d(0.0, theHeight), gp_Dir2d(1.0, 0.0)); + + const occ::handle aC2d = new Geom2dAdaptor_Curve(aPLine, aFirst, aLast); + const occ::handle aSurf = new GeomAdaptor_Surface(aCyl); + const occ::handle aCoS = new Adaptor3d_CurveOnSurface(aC2d, aSurf); + + GeomLib_CheckCurveOnSurface aCheck; + aCheck.Init(aC3d, Precision::PConfusion()); + aCheck.Perform(aCoS); + theMaxDistance = aCheck.MaxDistance(); + return aCheck; +} +} // namespace + +TEST(GeomLib_CheckCurveOnSurfaceTest, AnalyticCoincident_ReportsZeroDeviation) +{ + double aMaxDist = RealLast(); + const auto aCheck = makeCircleOnCylinderCheck(2.0, 2.0, 5.0, aMaxDist); + + EXPECT_TRUE(aCheck.IsDone()); + EXPECT_LE(aMaxDist, Precision::Confusion()); +} + +TEST(GeomLib_CheckCurveOnSurfaceTest, AnalyticDeviating_ReportsActualDeviation) +{ + // The 3D circle is 0.05 smaller than the cylinder carrying the pcurve: the + // deviation is constant and equal to the radius difference. This guards the + // sampled flat-deviation shortcut against swallowing real deviations. + double aMaxDist = RealLast(); + const auto aCheck = makeCircleOnCylinderCheck(1.95, 2.0, 5.0, aMaxDist); + + EXPECT_TRUE(aCheck.IsDone()); + EXPECT_NEAR(aMaxDist, 0.05, 1.0e-9); +} + +TEST(GeomLib_CheckCurveOnSurfaceTest, LineOnPlane_ReportsZeroDeviation) +{ + const double aFirst = -10.0, aLast = 10.0; + + const occ::handle aLine = new Geom_Line(gp_Pnt(0.0, 0.0, 0.0), gp_Dir(1.0, 0.0, 0.0)); + const occ::handle aC3d = new GeomAdaptor_Curve(aLine, aFirst, aLast); + + const occ::handle aPlane = + new Geom_Plane(gp_Ax3(gp_Pnt(0.0, 0.0, 0.0), gp_Dir(0.0, 0.0, 1.0), gp_Dir(1.0, 0.0, 0.0))); + const occ::handle aPLine = new Geom2d_Line(gp_Pnt2d(0.0, 0.0), gp_Dir2d(1.0, 0.0)); + + const occ::handle aC2d = new Geom2dAdaptor_Curve(aPLine, aFirst, aLast); + const occ::handle aSurf = new GeomAdaptor_Surface(aPlane); + const occ::handle aCoS = new Adaptor3d_CurveOnSurface(aC2d, aSurf); + + GeomLib_CheckCurveOnSurface aCheck; + aCheck.Init(aC3d, Precision::PConfusion()); + aCheck.Perform(aCoS); + + EXPECT_TRUE(aCheck.IsDone()); + EXPECT_LE(aCheck.MaxDistance(), Precision::Confusion()); +} diff --git a/src/ModelingData/TKGeomBase/GeomLib/GeomLib_CheckCurveOnSurface.cxx b/src/ModelingData/TKGeomBase/GeomLib/GeomLib_CheckCurveOnSurface.cxx index 695a0a409f..1bcf79c066 100644 --- a/src/ModelingData/TKGeomBase/GeomLib/GeomLib_CheckCurveOnSurface.cxx +++ b/src/ModelingData/TKGeomBase/GeomLib/GeomLib_CheckCurveOnSurface.cxx @@ -37,6 +37,7 @@ class GeomLib_CheckCurveOnSurface_TargetFunc; static bool MinComputing(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, const double theEpsilon, // 1.0e-3 const int theNbParticles, + const bool theIsAnalytic, double& theBestValue, double& theBestParameter); @@ -214,12 +215,14 @@ class GeomLib_CheckCurveOnSurface_Local const Array1OfHCurve& theCurveOnSurfaceArray, const NCollection_Array1& theIntervalsArr, const double theEpsilonRange, - const int theNbParticles) + const int theNbParticles, + const bool theIsAnalytic) : myCurveArray(theCurveArray), myCurveOnSurfaceArray(theCurveOnSurfaceArray), mySubIntervals(theIntervalsArr), myEpsilonRange(theEpsilonRange), myNbParticles(theNbParticles), + myIsAnalytic(theIsAnalytic), myArrOfDist(theIntervalsArr.Lower(), theIntervalsArr.Upper() - 1), myArrOfParam(theIntervalsArr.Lower(), theIntervalsArr.Upper() - 1) { @@ -239,7 +242,7 @@ class GeomLib_CheckCurveOnSurface_Local mySubIntervals.Value(theElemIndex + 1)); double aMinDist = RealLast(), aPar = 0.0; - if (!MinComputing(aFunc, myEpsilonRange, myNbParticles, aMinDist, aPar)) + if (!MinComputing(aFunc, myEpsilonRange, myNbParticles, myIsAnalytic, aMinDist, aPar)) { myArrOfDist(theElemIndex) = RealLast(); myArrOfParam(theElemIndex) = aFunc.FirstParameter(); @@ -278,6 +281,7 @@ class GeomLib_CheckCurveOnSurface_Local const NCollection_Array1& mySubIntervals; const double myEpsilonRange; const int myNbParticles; + const bool myIsAnalytic; mutable NCollection_Array1 myArrOfDist; mutable NCollection_Array1 myArrOfParam; }; @@ -332,6 +336,43 @@ void GeomLib_CheckCurveOnSurface::Init(const occ::handle& theCu //================================================================================================= +//================================================================================================= + +// Returns true when the whole composition is made of analytic geometry: the +// deviation function is then a low-degree smooth function whose behavior is +// fully resolved by the control-point sampling performed before the swarm +// optimization, so a flat sampled deviation can be trusted without running it. +static bool IsAnalyticComposition(const Adaptor3d_Curve& theCurve, + const Adaptor3d_CurveOnSurface& theCurveOnSurface) +{ + auto isAnalyticCurveType = [](const GeomAbs_CurveType theType) { + return theType == GeomAbs_Line || theType == GeomAbs_Circle || theType == GeomAbs_Ellipse + || theType == GeomAbs_Hyperbola || theType == GeomAbs_Parabola; + }; + + if (!isAnalyticCurveType(theCurve.GetType())) + { + return false; + } + + const occ::handle& aCurve2d = theCurveOnSurface.GetCurve(); + if (aCurve2d.IsNull() || !isAnalyticCurveType(aCurve2d->GetType())) + { + return false; + } + + const occ::handle& aSurface = theCurveOnSurface.GetSurface(); + if (aSurface.IsNull()) + { + return false; + } + const GeomAbs_SurfaceType aSurfType = aSurface->GetType(); + return aSurfType == GeomAbs_Plane || aSurfType == GeomAbs_Cylinder || aSurfType == GeomAbs_Cone + || aSurfType == GeomAbs_Sphere || aSurfType == GeomAbs_Torus; +} + +//================================================================================================= + void GeomLib_CheckCurveOnSurface::Perform( const occ::handle& theCurveOnSurface) { @@ -397,11 +438,14 @@ void GeomLib_CheckCurveOnSurface::Perform( aNbThreads > 1 ? theCurveOnSurface->ShallowCopy() : static_cast&>(theCurveOnSurface)); } + const bool isAnalytic = IsAnalyticComposition(*myCurve, *theCurveOnSurface); + GeomLib_CheckCurveOnSurface_Local aComp(aCurveArray, aCurveOnSurfaceArray, anIntervals, anEpsilonRange, - aNbParticles); + aNbParticles, + isAnalytic); if (aNbThreads > 1) { const occ::handle& aThreadPool = OSD_ThreadPool::DefaultPool(); @@ -646,9 +690,13 @@ bool PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, const math_Vector& theParSup, const double theEpsilon, const int theNbParticles, + const bool theIsAnalytic, double& theBestValue, - math_Vector& theOutputParam) + math_Vector& theOutputParam, + bool& theIsFlat) { + theIsFlat = false; + const double aDeltaParam = theParSup(1) - theParInf(1); if (aDeltaParam < Precision::PConfusion()) { @@ -663,6 +711,10 @@ bool PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, // They are used for finding a position of theNbParticles worst places const int aNbControlPoints = 3 * theNbParticles; + int aNbComputed = 0; + double aBestSeedVal = RealLast(); + double aBestSeedPrm = theParInf(1); + const double aStep = aDeltaParam / (aNbControlPoints - 1); int aCount = 1; for (double aPrm = theParInf(1); aCount <= aNbControlPoints; @@ -674,6 +726,13 @@ bool PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, continue; } + ++aNbComputed; + if (aVal < aBestSeedVal) + { + aBestSeedVal = aVal; + aBestSeedPrm = aPrm; + } + PSO_Particle* aParticle = aParticles.GetWorstParticle(); if (aVal > aParticle->BestDistance) @@ -687,6 +746,45 @@ bool PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, aParticle->BestDistance = aVal; } + // When every control point lies below the geometric noise floor (squared + // distances under Precision::SquareConfusion()), the curves are coincident + // within Precision::Confusion() and the swarm optimization, the Newton + // refinement (whose Hessian is singular on such a flat function) and the + // fallback swarm would only rediscover this value. For a fully analytic + // composition the control-point sampling already resolves the low-degree + // deviation function; otherwise flatness is confirmed by resampling at + // half-step offsets before the optimization is skipped. + if (aNbComputed == aNbControlPoints && aBestSeedVal >= -Precision::SquareConfusion()) + { + bool isConfirmedFlat = theIsAnalytic; + if (!isConfirmedFlat) + { + isConfirmedFlat = true; + for (double aPrm = theParInf(1) + 0.5 * aStep; aPrm < theParSup(1); aPrm += aStep) + { + double aVal = RealLast(); + if (!theFunction.Value(aPrm, aVal) || aVal < -Precision::SquareConfusion()) + { + isConfirmedFlat = false; + break; + } + if (aVal < aBestSeedVal) + { + aBestSeedVal = aVal; + aBestSeedPrm = aPrm; + } + } + } + + if (isConfirmedFlat) + { + theBestValue = aBestSeedVal; + theOutputParam(1) = aBestSeedPrm; + theIsFlat = true; + return true; + } + } + math_PSO aPSO(&theFunction, theParInf, theParSup, aStepPar); aPSO.Perform(aParticles, theNbParticles, theBestValue, theOutputParam); @@ -698,6 +796,7 @@ bool PSO_Perform(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, bool MinComputing(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, const double theEpsilon, // 1.0e-3 const int theNbParticles, + const bool theIsAnalytic, double& theBestValue, double& theBestParameter) { @@ -712,13 +811,16 @@ bool MinComputing(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, theBestParameter = aParInf(1); theBestValue = RealLast(); + bool isFlat = false; if (!PSO_Perform(theFunction, aParInf, aParSup, theEpsilon, theNbParticles, + theIsAnalytic, theBestValue, - anOutputParam)) + anOutputParam, + isFlat)) { #ifdef OCCT_DEBUG std::cout << "BRepLib_CheckCurveOnSurface::Compute(): math_PSO is failed!" << std::endl; @@ -728,6 +830,11 @@ bool MinComputing(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, theBestParameter = anOutputParam(1); + if (isFlat) + { + return true; + } + // Here, anOutputParam contains parameter, which is near to optimal. // It needs to be more precise. Precision is made by math_NewtonMinimum. math_NewtonMinimum aMinSol(theFunction); @@ -751,8 +858,10 @@ bool MinComputing(GeomLib_CheckCurveOnSurface_TargetFunc& theFunction, aParSup, theEpsilon, theNbParticles, + theIsAnalytic, aValue, - anOutputParam)) + anOutputParam, + isFlat)) { if (aValue < theBestValue) {