From bea50c2fd8b2b7f9e77fd59f38b7bebdc0557cd4 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Tue, 28 Apr 2026 12:50:12 +0100 Subject: [PATCH] Modeling - Tangent case with 2 cylinders (#1228) 2 cylinders are tangent along a line. The analytical intersection returns a point which is not valid and it is a restriction of the analytical intersection algorithm. The patch is ensured to use the implicit-implicit intersection algorithm. --- .../TKBO/GTests/FILES.cmake | 1 + .../TKBO/GTests/IntTools_FaceFace_Test.cxx | 424 ++++++++++++++++++ .../TKBool/GTests/FILES.cmake | 1 - .../TKBool/GTests/IntTools_FaceFace_Test.cxx | 61 --- .../IntPatch/IntPatch_ImpImpIntersection.cxx | 245 ++++------ .../GTests/IntAna_IntQuadQuad_Test.cxx | 67 +++ 6 files changed, 585 insertions(+), 214 deletions(-) create mode 100644 src/ModelingAlgorithms/TKBO/GTests/IntTools_FaceFace_Test.cxx delete mode 100644 src/ModelingAlgorithms/TKBool/GTests/IntTools_FaceFace_Test.cxx diff --git a/src/ModelingAlgorithms/TKBO/GTests/FILES.cmake b/src/ModelingAlgorithms/TKBO/GTests/FILES.cmake index 9430a60054..e78fe1d017 100644 --- a/src/ModelingAlgorithms/TKBO/GTests/FILES.cmake +++ b/src/ModelingAlgorithms/TKBO/GTests/FILES.cmake @@ -9,4 +9,5 @@ set(OCCT_TKBO_GTests_FILES BRepAlgoAPI_Common_Test.cxx BOPAlgo_BOP_Test.cxx BOPAlgo_PaveFiller_Test.cxx + IntTools_FaceFace_Test.cxx ) diff --git a/src/ModelingAlgorithms/TKBO/GTests/IntTools_FaceFace_Test.cxx b/src/ModelingAlgorithms/TKBO/GTests/IntTools_FaceFace_Test.cxx new file mode 100644 index 0000000000..4489a5806c --- /dev/null +++ b/src/ModelingAlgorithms/TKBO/GTests/IntTools_FaceFace_Test.cxx @@ -0,0 +1,424 @@ +// Copyright (c) 2025 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 +#include +#include +#include +#include +#include +#include +#include +#include + +//================================================================================================== +// IntTools_FaceFace Tests - Tests for face-face intersection +//================================================================================================== + +class IntTools_FaceFaceTest : public ::testing::Test +{ +protected: + //! Create a half-cylinder face. + //! @param theCenter center point on the cylinder axis + //! @param theAxis direction of the cylinder axis + //! @param theRadius radius of the cylinder + //! @param theVMin minimum V parameter (height along axis) + //! @param theVMax maximum V parameter (height along axis) + //! @param theUMin minimum U parameter (angle, default 0) + //! @param theUMax maximum U parameter (angle, default PI for half-cylinder) + //! @return the half-cylinder face + static TopoDS_Face CreateHalfCylinderFace(const gp_Pnt& theCenter, + const gp_Dir& theAxis, + double theRadius, + double theVMin, + double theVMax, + double theUMin = 0.0, + double theUMax = M_PI) + { + // Create cylinder surface + gp_Ax3 anAx3(theCenter, theAxis); + Handle(Geom_CylindricalSurface) aCylSurf = new Geom_CylindricalSurface(anAx3, theRadius); + + // Create face with parameter bounds + BRepBuilderAPI_MakeFace aFaceMaker(aCylSurf, theUMin, theUMax, theVMin, theVMax, 1e-7); + return aFaceMaker.Face(); + } + + //! Create a full cylinder face with finite V range. + //! @param theCenter center point on the cylinder axis + //! @param theAxis direction of the cylinder axis + //! @param theRadius radius of the cylinder + //! @param theVMin minimum V parameter (height along axis) + //! @param theVMax maximum V parameter (height along axis) + //! @return the cylinder face + static TopoDS_Face CreateCylinderFace(const gp_Pnt& theCenter, + const gp_Dir& theAxis, + double theRadius, + double theVMin, + double theVMax) + { + return CreateHalfCylinderFace(theCenter, theAxis, theRadius, theVMin, theVMax, 0.0, 2.0 * M_PI); + } + + //! Create a circular planar face. + //! @param theCenter center of the circle on the plane + //! @param theNormal normal direction of the plane + //! @param theRadius radius of the circular boundary + //! @return the circular planar face + static TopoDS_Face CreateCircularPlaneFace(const gp_Pnt& theCenter, + const gp_Dir& theNormal, + double theRadius) + { + // Create plane + gp_Pln aPln(theCenter, theNormal); + Handle(Geom_Plane) aPlane = new Geom_Plane(aPln); + + // Create circular edge on the plane + gp_Ax2 aCircAx(theCenter, theNormal); + gp_Circ aCirc(aCircAx, theRadius); + + BRepBuilderAPI_MakeEdge anEdgeMaker(aCirc); + TopoDS_Edge aCircEdge = anEdgeMaker.Edge(); + + // Create wire from circular edge + BRepBuilderAPI_MakeWire aWireMaker(aCircEdge); + TopoDS_Wire aWire = aWireMaker.Wire(); + + // Create face from plane with circular wire + BRepBuilderAPI_MakeFace aFaceMaker(aPlane, aWire, true); + return aFaceMaker.Face(); + } +}; + +//! Regression: cylinder-cylinder analytical gating must not depend on argument order. +//! Geometry is configured so that only one cylinder contributes a touching V-boundary. +TEST_F(IntTools_FaceFaceTest, PerpendicularCylinderBoundaryTouch_OrderIndependent) +{ + // Cylinder 1: axis along Z, long finite segment. + const gp_Pnt aCyl1Center(0.0, 0.0, 0.0); + const gp_Dir aCyl1Axis(0.0, 0.0, 1.0); + const double aCyl1Radius = 2.0; + const double aCyl1VMin = -5.0; + const double aCyl1VMax = 5.0; + + // Cylinder 2: axis along X, one V-boundary at X=0 touches cylinder 1. + // Endpoint (0, 2, 0) lies on cylinder 1 surface, while cylinder 1 V-boundaries + // stay far from cylinder 2 surface. + const gp_Pnt aCyl2Center(0.0, 2.0, 0.0); + const gp_Dir aCyl2Axis(1.0, 0.0, 0.0); + const double aCyl2Radius = 0.5; + const double aCyl2VMin = 0.0; + const double aCyl2VMax = 4.0; + + TopoDS_Face aFace1 = + CreateCylinderFace(aCyl1Center, aCyl1Axis, aCyl1Radius, aCyl1VMin, aCyl1VMax); + TopoDS_Face aFace2 = + CreateCylinderFace(aCyl2Center, aCyl2Axis, aCyl2Radius, aCyl2VMin, aCyl2VMax); + + ASSERT_FALSE(aFace1.IsNull()) << "Failed to create first cylinder face"; + ASSERT_FALSE(aFace2.IsNull()) << "Failed to create second cylinder face"; + + IntTools_FaceFace aFF12; + aFF12.SetParameters(true, true, true, 1.0e-7); + aFF12.Perform(aFace1, aFace2); + ASSERT_TRUE(aFF12.IsDone()) << "Intersection failed for (face1, face2)"; + + IntTools_FaceFace aFF21; + aFF21.SetParameters(true, true, true, 1.0e-7); + aFF21.Perform(aFace2, aFace1); + ASSERT_TRUE(aFF21.IsDone()) << "Intersection failed for (face2, face1)"; + + const int aNbCurves12 = aFF12.Lines().Length(); + const int aNbPoints12 = aFF12.Points().Length(); + const int aNbCurves21 = aFF21.Lines().Length(); + const int aNbPoints21 = aFF21.Points().Length(); + + EXPECT_EQ(aNbCurves12, aNbCurves21) + << "Intersection curve count should not depend on argument order"; + EXPECT_EQ(aNbPoints12, aNbPoints21) + << "Intersection point count should not depend on argument order"; +} + +//! Test that a half-cylinder face and a circular plane face that don't actually +//! intersect (cylinder is outside the circle) correctly return no intersection. +//! This tests the fix for wire-based boundary validation in IntTools_FaceFace. +TEST_F(IntTools_FaceFaceTest, HalfCylinderOutsideCircularPlane_NoIntersection) +{ + // Create geometry similar to the reported bug: + // - Half cylinder at position where it would intersect the infinite plane + // - But the circular boundary of the plane doesn't reach the cylinder + + // Plane at Y = 0, with circular boundary of radius 10, centered at origin + const gp_Pnt aPlaneCenter(0.0, 0.0, 0.0); + const gp_Dir aPlaneNormal(0.0, 1.0, 0.0); + const double aPlaneRadius = 10.0; + + // Half cylinder with axis along Y, positioned far from plane center + // The cylinder is at X=20, which is outside the plane's circular boundary (radius 10) + const gp_Pnt aCylCenter(20.0, 0.0, 0.0); + const gp_Dir aCylAxis(0.0, 1.0, 0.0); + const double aCylRadius = 2.0; + const double aVMin = -1.0; + const double aVMax = 1.0; + + // Create faces + TopoDS_Face aCircularPlane = CreateCircularPlaneFace(aPlaneCenter, aPlaneNormal, aPlaneRadius); + TopoDS_Face aHalfCylinder = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax); + + ASSERT_FALSE(aCircularPlane.IsNull()) << "Failed to create circular plane face"; + ASSERT_FALSE(aHalfCylinder.IsNull()) << "Failed to create half-cylinder face"; + + // Perform face-face intersection + IntTools_FaceFace aFF; + aFF.SetParameters(true, true, true, 1.0e-7); + aFF.Perform(aCircularPlane, aHalfCylinder); + + ASSERT_TRUE(aFF.IsDone()) << "IntTools_FaceFace::Perform() failed"; + + // The intersection should be empty because the cylinder is outside the circular boundary + const int aNbCurves = aFF.Lines().Length(); + const int aNbPoints = aFF.Points().Length(); + + EXPECT_EQ(aNbCurves, 0) << "Expected no intersection curves, but found " << aNbCurves + << ". The cylinder at X=20 is outside the circular plane (radius 10)."; + EXPECT_EQ(aNbPoints, 0) << "Expected no intersection points, but found " << aNbPoints; +} + +//! Test that a half-cylinder face and a circular plane face that do intersect +//! correctly return intersection curves. +TEST_F(IntTools_FaceFaceTest, HalfCylinderInsideCircularPlane_HasIntersection) +{ + // Create geometry where intersection should exist: + // - Half cylinder positioned inside the circular plane boundary + + // Plane at Y = 0, with circular boundary of radius 20, centered at origin + const gp_Pnt aPlaneCenter(0.0, 0.0, 0.0); + const gp_Dir aPlaneNormal(0.0, 1.0, 0.0); + const double aPlaneRadius = 20.0; + + // Half cylinder with axis along Y, positioned at origin (inside plane circle) + const gp_Pnt aCylCenter(0.0, 0.0, 0.0); + const gp_Dir aCylAxis(0.0, 1.0, 0.0); + const double aCylRadius = 5.0; + const double aVMin = -1.0; + const double aVMax = 1.0; + + // Create faces + TopoDS_Face aCircularPlane = CreateCircularPlaneFace(aPlaneCenter, aPlaneNormal, aPlaneRadius); + TopoDS_Face aHalfCylinder = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax); + + ASSERT_FALSE(aCircularPlane.IsNull()) << "Failed to create circular plane face"; + ASSERT_FALSE(aHalfCylinder.IsNull()) << "Failed to create half-cylinder face"; + + // Perform face-face intersection + IntTools_FaceFace aFF; + aFF.SetParameters(true, true, true, 1.0e-7); + aFF.Perform(aCircularPlane, aHalfCylinder); + + ASSERT_TRUE(aFF.IsDone()) << "IntTools_FaceFace::Perform() failed"; + + // The intersection should exist - a semicircular arc + const int aNbCurves = aFF.Lines().Length(); + + EXPECT_GE(aNbCurves, 1) + << "Expected at least one intersection curve (semicircular arc), but found " << aNbCurves; +} + +//! Test geometry where the opposite half of the cylinder IS inside the circle. +//! Position the cylinder closer so that the half bulging toward the plane is inside the circle. +TEST_F(IntTools_FaceFaceTest, OppositeHalfInsideCircle_HasIntersection) +{ + // Plane circle centered at origin with radius 15 + const gp_Pnt aPlaneCenter(0.0, 0.0, 0.0); + const gp_Dir aPlaneNormal(0.0, 1.0, 0.0); + const double aPlaneRadius = 15.0; + + // Cylinder at X=10, Z=-5 + // The half with U in [PI, 2*PI] bulges toward +Z (toward origin) + // At U=3*PI/2, the point is at (10, Y, -5 + 3) = (10, Y, -2) + // Distance from (10, -2) to origin = sqrt(100 + 4) ~ 10.2 < 15 (inside!) + const gp_Pnt aCylCenter(10.0, 0.0, -5.0); + const gp_Dir aCylAxis(0.0, 1.0, 0.0); + const double aCylRadius = 3.0; + const double aVMin = -1.0; + const double aVMax = 1.0; + + // Create the opposite half: U in [PI, 2*PI] (bulges toward +Z, toward plane center) + TopoDS_Face aCircularPlane = CreateCircularPlaneFace(aPlaneCenter, aPlaneNormal, aPlaneRadius); + TopoDS_Face aHalfCylinder = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax, M_PI, 2.0 * M_PI); + + ASSERT_FALSE(aCircularPlane.IsNull()) << "Failed to create circular plane face"; + ASSERT_FALSE(aHalfCylinder.IsNull()) << "Failed to create half-cylinder face"; + + // Perform face-face intersection + IntTools_FaceFace aFF; + aFF.SetParameters(true, true, true, 1.0e-7); + aFF.Perform(aCircularPlane, aHalfCylinder); + + ASSERT_TRUE(aFF.IsDone()) << "IntTools_FaceFace::Perform() failed"; + + // This intersection SHOULD exist - the semicircle is partially inside the plane circle + const int aNbCurves = aFF.Lines().Length(); + + EXPECT_GE(aNbCurves, 1) + << "Expected intersection for opposite half that bulges toward plane center."; +} + +//! Test geometry where the cylinder arc partially crosses the circle boundary. +//! With wire-aware BRepTopAdaptor_TopolTool, the intersection curve is properly +//! trimmed to the portion that lies within both face wire boundaries. +TEST_F(IntTools_FaceFaceTest, PartialCrossing_ProperlyTrimmed) +{ + const gp_Pnt aPlaneCenter(0.0, 0.0, 0.0); + const gp_Dir aPlaneNormal(0.0, 1.0, 0.0); + const double aPlaneRadius = 15.0; + + // Position cylinder so the arc partially crosses the boundary. + // First half: U in [0, PI], bulges toward -Z + // At U=PI/2: (12, Y, -11), distance = sqrt(144 + 121) ~ 16.3 > 15 (outside!) + // At U=0: (15, Y, -8), distance = sqrt(225 + 64) ~ 17.0 > 15 (outside!) + // At U=PI: (9, Y, -8), distance = sqrt(81 + 64) ~ 12.0 < 15 (inside!) + // The arc crosses from outside to inside the circle. + const gp_Pnt aCylCenter(12.0, 0.0, -8.0); + const gp_Dir aCylAxis(0.0, 1.0, 0.0); + const double aCylRadius = 3.0; + const double aVMin = -1.0; + const double aVMax = 1.0; + + TopoDS_Face aCircularPlane = CreateCircularPlaneFace(aPlaneCenter, aPlaneNormal, aPlaneRadius); + TopoDS_Face aHalfCylinder = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax, 0.0, M_PI); + + ASSERT_FALSE(aCircularPlane.IsNull()) << "Failed to create circular plane face"; + ASSERT_FALSE(aHalfCylinder.IsNull()) << "Failed to create half-cylinder face"; + + IntTools_FaceFace aFF; + aFF.SetParameters(true, true, true, 1.0e-7); + aFF.Perform(aCircularPlane, aHalfCylinder); + + ASSERT_TRUE(aFF.IsDone()) << "IntTools_FaceFace::Perform() failed"; + + const int aNbCurves = aFF.Lines().Length(); + + // With wire-aware BRepTopAdaptor_TopolTool, partial intersections are properly + // trimmed to the valid region within both face boundaries, not rejected entirely. + // The portion of the cylinder that's inside the circular plane produces a valid + // intersection curve segment. + EXPECT_GE(aNbCurves, 1) + << "Wire-aware classification should produce trimmed intersection curve for partial overlap"; +} + +//! Test with cylinder positioned such that both halves are completely outside the circle. +//! This verifies the algorithm correctly handles cases where no part of the cylinder +//! is inside the circular boundary. +TEST_F(IntTools_FaceFaceTest, BothHalvesCompletelyOutside_NoIntersection) +{ + // Plane circle with small radius + const gp_Pnt aPlaneCenter(0.0, 0.0, 0.0); + const gp_Dir aPlaneNormal(0.0, 1.0, 0.0); + const double aPlaneRadius = 5.0; + + // Cylinder far from the plane center + // Both halves will be completely outside the circle + const gp_Pnt aCylCenter(20.0, 0.0, 20.0); + const gp_Dir aCylAxis(0.0, 1.0, 0.0); + const double aCylRadius = 2.0; + const double aVMin = -1.0; + const double aVMax = 1.0; + + // Test first half + TopoDS_Face aCircularPlane = CreateCircularPlaneFace(aPlaneCenter, aPlaneNormal, aPlaneRadius); + TopoDS_Face aHalfCylinder1 = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax, 0.0, M_PI); + + ASSERT_FALSE(aCircularPlane.IsNull()); + ASSERT_FALSE(aHalfCylinder1.IsNull()); + + IntTools_FaceFace aFF1; + aFF1.SetParameters(true, true, true, 1.0e-7); + aFF1.Perform(aCircularPlane, aHalfCylinder1); + + ASSERT_TRUE(aFF1.IsDone()); + EXPECT_EQ(aFF1.Lines().Length(), 0) << "First half should have no intersection"; + EXPECT_EQ(aFF1.Points().Length(), 0); + + // Test opposite half + TopoDS_Face aHalfCylinder2 = + CreateHalfCylinderFace(aCylCenter, aCylAxis, aCylRadius, aVMin, aVMax, M_PI, 2.0 * M_PI); + + ASSERT_FALSE(aHalfCylinder2.IsNull()); + + IntTools_FaceFace aFF2; + aFF2.SetParameters(true, true, true, 1.0e-7); + aFF2.Perform(aCircularPlane, aHalfCylinder2); + + ASSERT_TRUE(aFF2.IsDone()); + EXPECT_EQ(aFF2.Lines().Length(), 0) << "Opposite half should also have no intersection"; + EXPECT_EQ(aFF2.Points().Length(), 0); +} + +// Test OCC24005: IntTools_FaceFace should complete quickly and correctly when +// intersecting a slightly off-angle plane with a cylinder. Previously took 7+ seconds. +TEST_F(IntTools_FaceFaceTest, OCC24005_PlaneCylinderIntersection) +{ + // Hardcoded geometry from the original regression test + Handle(Geom_Plane) aPlane = new Geom_Plane( + gp_Ax3(gp_Pnt(-72.948737453424499, 754.30437716359393, 259.52151854671678), + gp_Dir(6.2471473085930200e-007, -0.99999999999980493, 0.00000000000000000), + gp_Dir(0.99999999999980493, 6.2471473085930200e-007, 0.00000000000000000))); + + Handle(Geom_CylindricalSurface) aCylinder = new Geom_CylindricalSurface( + gp_Ax3(gp_Pnt(-6.4812490053250649, 753.39408794522092, 279.16400974257465), + gp_Dir(1.0000000000000000, 0.0, 0.00000000000000000), + gp_Dir(0.0, 1.0000000000000000, 0.00000000000000000)), + 19.712534607908712); + + BRep_Builder aBuilder; + TopoDS_Face aFace1, aFace2; + aBuilder.MakeFace(aFace1, aPlane, Precision::Confusion()); + aBuilder.MakeFace(aFace2, aCylinder, Precision::Confusion()); + + IntTools_FaceFace anInters; + anInters.SetParameters(false, true, true, Precision::Confusion()); + anInters.Perform(aFace1, aFace2); + + ASSERT_TRUE(anInters.IsDone()) << "IntTools_FaceFace::Perform did not complete"; + + // The original test verified that intersection completes without hanging. + // Check that we get at least one intersection curve. + const NCollection_Sequence& aCurves = anInters.Lines(); + EXPECT_GE(aCurves.Length() + anInters.Points().Length(), 1) + << "Expected at least one intersection result (curve or point)"; +} diff --git a/src/ModelingAlgorithms/TKBool/GTests/FILES.cmake b/src/ModelingAlgorithms/TKBool/GTests/FILES.cmake index 35beaf4fdd..9e9735b43c 100644 --- a/src/ModelingAlgorithms/TKBool/GTests/FILES.cmake +++ b/src/ModelingAlgorithms/TKBool/GTests/FILES.cmake @@ -6,5 +6,4 @@ set(OCCT_TKBool_GTests_FILES BRepAlgoAPI_Fuse_Test.cxx BRepAlgoAPI_Section_Test.cxx BRepFill_PipeShell_Test.cxx - IntTools_FaceFace_Test.cxx ) diff --git a/src/ModelingAlgorithms/TKBool/GTests/IntTools_FaceFace_Test.cxx b/src/ModelingAlgorithms/TKBool/GTests/IntTools_FaceFace_Test.cxx deleted file mode 100644 index ba5e3ff62d..0000000000 --- a/src/ModelingAlgorithms/TKBool/GTests/IntTools_FaceFace_Test.cxx +++ /dev/null @@ -1,61 +0,0 @@ -// 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 - -// Test OCC24005: IntTools_FaceFace should complete quickly and correctly when -// intersecting a slightly off-angle plane with a cylinder. Previously took 7+ seconds. -TEST(IntTools_FaceFaceTest, OCC24005_PlaneCylinderIntersection) -{ - // Hardcoded geometry from the original regression test - Handle(Geom_Plane) aPlane = new Geom_Plane( - gp_Ax3(gp_Pnt(-72.948737453424499, 754.30437716359393, 259.52151854671678), - gp_Dir(6.2471473085930200e-007, -0.99999999999980493, 0.00000000000000000), - gp_Dir(0.99999999999980493, 6.2471473085930200e-007, 0.00000000000000000))); - - Handle(Geom_CylindricalSurface) aCylinder = new Geom_CylindricalSurface( - gp_Ax3(gp_Pnt(-6.4812490053250649, 753.39408794522092, 279.16400974257465), - gp_Dir(1.0000000000000000, 0.0, 0.00000000000000000), - gp_Dir(0.0, 1.0000000000000000, 0.00000000000000000)), - 19.712534607908712); - - BRep_Builder aBuilder; - TopoDS_Face aFace1, aFace2; - aBuilder.MakeFace(aFace1, aPlane, Precision::Confusion()); - aBuilder.MakeFace(aFace2, aCylinder, Precision::Confusion()); - - IntTools_FaceFace anInters; - anInters.SetParameters(false, true, true, Precision::Confusion()); - anInters.Perform(aFace1, aFace2); - - ASSERT_TRUE(anInters.IsDone()) << "IntTools_FaceFace::Perform did not complete"; - - // The original test verified that intersection completes without hanging. - // Check that we get at least one intersection curve. - const NCollection_Sequence& aCurves = anInters.Lines(); - EXPECT_GE(aCurves.Length() + anInters.Points().Length(), 1) - << "Expected at least one intersection result (curve or point)"; -} diff --git a/src/ModelingAlgorithms/TKGeomAlgo/IntPatch/IntPatch_ImpImpIntersection.cxx b/src/ModelingAlgorithms/TKGeomAlgo/IntPatch/IntPatch_ImpImpIntersection.cxx index f94683bf32..05c45bf9fc 100644 --- a/src/ModelingAlgorithms/TKGeomAlgo/IntPatch/IntPatch_ImpImpIntersection.cxx +++ b/src/ModelingAlgorithms/TKGeomAlgo/IntPatch/IntPatch_ImpImpIntersection.cxx @@ -51,6 +51,7 @@ #include #include #include +#include #include #include @@ -4210,6 +4211,10 @@ ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1, mC /= aA; } +// Cyl-cyl-specific helper: constructed only by IntCyCy (see ~line 8495/8517) +// and relies on cylinder-cylinder coefficient tables (ComputationMethods:: +// stCoeffsValue) for both SearchOnVBounds' Newton system and the analytical +// CylCylComputeParameters evaluators used downstream. class WorkWithBoundaries { public: @@ -4295,8 +4300,13 @@ public: // Returns UV-bounds of 2nd surface const Bnd_Box2d& UVS2() const { return myUVSurf2; } + // theU1Prev is the walker's previous U1 sample. Together with theU1 it + // forms the bracket across which the straddle on V_bound was detected and + // in which the V-boundary crossing is resolved by 1D branch-aware + // bisection on the analytical intersection curve. void AddBoundaryPoint(const occ::handle& theWL, const double theU1, + const double theU1Prev, const double theU1Min, const double theU2, const double theV1, @@ -4318,17 +4328,6 @@ public: Bnd_Range& theOutBoxS2) const; protected: - // Solves equation (2) (see declaration of ComputationMethods class) in case, - // when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary - // and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and - // requires initial values (theVInit, theInitU2 and theInitMainVar). - bool SearchOnVBounds(const SearchBoundType theSBType, - const double theVzad, - const double theVInit, - const double theInitU2, - const double theInitMainVar, - double& theMainVariableValue) const; - const WorkWithBoundaries& operator=(const WorkWithBoundaries&); private: @@ -5386,129 +5385,6 @@ bool ComputationMethods::CylCylComputeParameters(const double theU1par, return true; } -//================================================================================================= - -bool WorkWithBoundaries::SearchOnVBounds(const SearchBoundType theSBType, - const double theVzad, - const double theVInit, - const double theInitU2, - const double theInitMainVar, - double& theMainVariableValue) const -{ - const int aNbDim = 3; - const double aMaxError = 4.0 * M_PI; // two periods - - theMainVariableValue = theInitMainVar; - const double aTol2 = 1.0e-18; - double aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit; - - // Structure of aMatr: - // C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*}, - // where C_{1}, C_{2} and C_{3} are math_Vector. - math_Matrix aMatr(1, aNbDim, 1, aNbDim); - - double anError = RealLast(); - double anErrorPrev = anError; - int aNbIter = 0; - do - { - if (++aNbIter > 1000) - return false; - - const double aSinU1 = sin(aMainVarPrev), aCosU1 = cos(aMainVarPrev), aSinU2 = sin(aU2Prev), - aCosU2 = cos(aU2Prev); - - math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev + myCoeffs.mVecB2) * aSinU2 - - (myCoeffs.mVecB2 * aU2Prev - myCoeffs.mVecA2) * aCosU2 - + (myCoeffs.mVecA1 * aMainVarPrev + myCoeffs.mVecB1) * aSinU1 - - (myCoeffs.mVecB1 * aMainVarPrev - myCoeffs.mVecA1) * aCosU1 - + myCoeffs.mVecD; - - math_Vector aMSum(1, 3); - - switch (theSBType) - { - case SearchV1: - aMatr.SetCol(3, myCoeffs.mVecC2); - aMSum = myCoeffs.mVecC1 * theVzad; - aVecFreeMem -= aMSum; - aMSum += myCoeffs.mVecC2 * anOtherVar; - break; - - case SearchV2: - aMatr.SetCol(3, myCoeffs.mVecC1); - aMSum = myCoeffs.mVecC2 * theVzad; - aVecFreeMem -= aMSum; - aMSum += myCoeffs.mVecC1 * anOtherVar; - break; - - default: - return false; - } - - aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1); - aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2); - - double aDetMainSyst = aMatr.Determinant(); - - if (std::abs(aDetMainSyst) < aNulValue) - { - return false; - } - - math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr); - aM1.SetCol(1, aVecFreeMem); - aM2.SetCol(2, aVecFreeMem); - aM3.SetCol(3, aVecFreeMem); - - const double aDetMainVar = aM1.Determinant(); - const double aDetVar1 = aM2.Determinant(); - const double aDetVar2 = aM3.Determinant(); - - double aDelta = aDetMainVar / aDetMainSyst - aMainVarPrev; - - if (std::abs(aDelta) > aMaxError) - return false; - - anError = aDelta * aDelta; - aMainVarPrev += aDelta; - - /// - aDelta = aDetVar1 / aDetMainSyst - aU2Prev; - - if (std::abs(aDelta) > aMaxError) - return false; - - anError += aDelta * aDelta; - aU2Prev += aDelta; - - /// - aDelta = aDetVar2 / aDetMainSyst - anOtherVar; - anError += aDelta * aDelta; - anOtherVar += aDelta; - - if (anError > anErrorPrev) - { // Method diverges. Keep the best result - const double aSinU1Last = sin(aMainVarPrev), aCosU1Last = cos(aMainVarPrev), - aSinU2Last = sin(aU2Prev), aCosU2Last = cos(aU2Prev); - aMSum -= (myCoeffs.mVecA1 * aCosU1Last + myCoeffs.mVecB1 * aSinU1Last - + myCoeffs.mVecA2 * aCosU2Last + myCoeffs.mVecB2 * aSinU2Last + myCoeffs.mVecD); - const double aSQNorm = aMSum.Norm2(); - return (aSQNorm < aTol2); - } - else - { - theMainVariableValue = aMainVarPrev; - } - - anErrorPrev = anError; - } while (anError > aTol2); - - theMainVariableValue = aMainVarPrev; - - return true; -} - //======================================================================= // function : InscribePoint // purpose : If theFlForce==TRUE theUGiven will be changed forcefully @@ -5859,6 +5735,7 @@ static bool AddPointIntoWL(const IntSurf_Quadric& theQuad1, //======================================================================= void WorkWithBoundaries::AddBoundaryPoint(const occ::handle& theWL, const double theU1, + const double theU1Prev, const double theU1Min, const double theU2, const double theV1, @@ -5883,6 +5760,65 @@ void WorkWithBoundaries::AddBoundaryPoint(const occ::handle& the StPInfo aUVPoint[aSize]; + // Branch-aware 1D root-finder for V(U1) = V_bound on the analytical + // intersection curve, used in place of the original rank-deficient 3x3 + // Newton (SearchOnVBounds). The intersection curve is 1D in U1 on each + // arccos branch, so the correct numerical tool is a bracketed 1D root + // finder. We use MathRoot::Brent, which converges robustly even when the + // curve grazes V_bound at a near-tangent (two close crossings on opposite + // sides of a V-extremum) -- the configuration that makes the 3x3 Newton + // Jacobian rank-deficient. The only tolerance is Precision::PConfusion() + // (OCCT standard parametric equality) for termination on U. + struct VBoundDelta + { + const ComputationMethods::stCoeffsValue& Coeffs; + int WLIndex; + bool IsV1; + double VBound; + + bool Value(double theX, double& theF) const + { + double aU2 = 0.0, aV1 = 0.0, aV2 = 0.0; + if (!ComputationMethods::CylCylComputeParameters(theX, WLIndex, Coeffs, aU2, aV1, aV2)) + return false; + theF = (IsV1 ? aV1 : aV2) - VBound; + return true; + } + }; + + auto findVBoundCrossing = [this, theWLIndex](const bool theIsV1, + const double theVBound, + const double theULo, + const double theUHi, + double& theU1Star) -> bool { + if (!(theULo < theUHi)) + return false; + VBoundDelta aFunc{myCoeffs, theWLIndex, theIsV1, theVBound}; + double aFlo = 0.0, aFhi = 0.0; + if (!aFunc.Value(theULo, aFlo) || !aFunc.Value(theUHi, aFhi)) + return false; + if (aFlo == 0.0) + { + theU1Star = theULo; + return true; + } + if (aFhi == 0.0) + { + theU1Star = theUHi; + return true; + } + if (aFlo * aFhi > 0.0) + return false; // no bracketed crossing + + MathUtils::Config aCfg; + aCfg.XTolerance = Precision::PConfusion(); + const MathUtils::ScalarResult aRes = MathRoot::Brent(aFunc, theULo, theUHi, aCfg); + if (aRes.Status != MathUtils::Status::OK || !aRes.Root.has_value()) + return false; + theU1Star = *aRes.Root; + return true; + }; + for (int anIDSurf = 0; anIDSurf < 4; anIDSurf += 2) { const double aVf = (anIDSurf == 0) ? theV1 : theV2, @@ -5902,16 +5838,12 @@ void WorkWithBoundaries::AddBoundaryPoint(const occ::handle& the continue; } - // Segment [aVf, aVl] intersects at least one V-boundary (first or last) - // (in general, case is possible, when aVf > aVl). - - // Precise intersection point - const bool aRes = SearchOnVBounds(aTS, - anArrVzad[anIndex], - (anIDSurf == 0) ? theV2 : theV1, - theU2, - theU1, - aUVPoint[anIndex].myU1); + // Segment [aVf, aVl] intersects at least one V-boundary (first or last). + // Resolve the crossing by 1D bisection on the analytical curve. + const double aULo = std::min(theU1Prev, theU1); + const double aUHi = std::max(theU1Prev, theU1); + const bool aRes = + findVBoundCrossing(aTS == SearchV1, anArrVzad[anIndex], aULo, aUHi, aUVPoint[anIndex].myU1); // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than // to theU1+/-Period @@ -6432,8 +6364,6 @@ static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& the // Out of "while(theNbCritPointsMax > 0)" cycle. break; } - - // Attention! Here theU1crit may be unsorted. } //======================================================================= @@ -6657,10 +6587,9 @@ static IntPatch_ImpImpIntersection::IntStatus CyCyNoGeometric( Precision::Infinite(), Precision::Infinite()}; - // This list of critical points is not full because it does not contain any points - // which intersection line goes through V-bounds of cylinders in. - // They are computed by numerical methods on - line (during algorithm working). - // The moment is caught, when intersection line goes through V-bounds of any cylinder. + // V-boundary points (both crossings and tangencies) are caught on-the-fly + // by the walker via AddBoundaryPoint; see the 1D branch-aware root-finder + // it uses in place of the rank-deficient 3x3 Newton SearchOnVBounds. int aNbCritPoints = aNbCritPointsMax; CriticalPointsComputing(anEquationCoeffs, @@ -6743,7 +6672,8 @@ static IntPatch_ImpImpIntersection::IntStatus CyCyNoGeometric( } double anU1 = anUf, aMinCriticalParam = anUf; - bool isFirst = true; + double anU1Prev = anUf; + bool isFirst = true; while (anU1 <= anUl) { @@ -6974,6 +6904,7 @@ static IntPatch_ImpImpIntersection::IntStatus CyCyNoGeometric( theBW.AddBoundaryPoint(aWLine[i], anU1, + anU1Prev, aMinCriticalParam, aU2[i], aV1[i], @@ -7098,6 +7029,7 @@ static IntPatch_ImpImpIntersection::IntStatus CyCyNoGeometric( theBW.AddBoundaryPoint(aWLine[i], anU1, + anU1Prev, aMinCriticalParam, aU2[i], aV1[i], @@ -7326,7 +7258,8 @@ static IntPatch_ImpImpIntersection::IntStatus CyCyNoGeometric( aMinUexp = std::min(aMinUexp, anUexpect[i]); } - anU1 = aMinUexp; + anU1Prev = anU1; + anU1 = aMinUexp; } if (Precision::PConfusion() >= (anUl - anU1)) @@ -7779,6 +7712,14 @@ IntPatch_ImpImpIntersection::IntStatus IntCyCy( { return IntPatch_ImpImpIntersection::IntStatus_OK; } + + // Analytical handler declined: discard anything it may have appended + // or flipped so the numerical path starts from a clean output state. + isTheEmpty = true; + isTheSameSurface = false; + isTheMultiplePoint = false; + theSlin.Clear(); + theSPnt.Clear(); } // Here, intersection line is not an analytical curve(line, circle, ellipsis etc.) diff --git a/src/ModelingData/TKGeomBase/GTests/IntAna_IntQuadQuad_Test.cxx b/src/ModelingData/TKGeomBase/GTests/IntAna_IntQuadQuad_Test.cxx index ead0184079..6ed2a9dc59 100644 --- a/src/ModelingData/TKGeomBase/GTests/IntAna_IntQuadQuad_Test.cxx +++ b/src/ModelingData/TKGeomBase/GTests/IntAna_IntQuadQuad_Test.cxx @@ -12,12 +12,14 @@ // commercial license or contractual agreement. #include +#include #include #include #include #include #include #include +#include #include @@ -204,3 +206,68 @@ TEST_F(IntAna_IntQuadQuad_Test, IndexingConsistencyTest) } } } + +// Regression test for near-tangent parallel cylinders. +// Very small chord between two theoretical lines must collapse to one line. +TEST_F(IntAna_IntQuadQuad_Test, CylinderCylinderNearTangent_CollapsesToSingleLine) +{ + gp_Pnt aCenter1(0.0, 0.0, 0.0); + gp_Pnt aCenter2(2.0 - 1.0e-9, 0.0, 0.0); + gp_Dir aZDir(gp_Dir::D::Z); + gp_Dir anXDir(gp_Dir::D::X); + gp_Ax3 anAxis1(aCenter1, aZDir, anXDir); + gp_Ax3 anAxis2(aCenter2, aZDir, anXDir); + + gp_Cylinder aCyl1(anAxis1, 1.0); + gp_Cylinder aCyl2(anAxis2, 1.0); + + IntAna_QuadQuadGeo anInter(aCyl1, aCyl2, 1.0e-4); + + EXPECT_TRUE(anInter.IsDone()); + EXPECT_EQ(anInter.TypeInter(), IntAna_Line); + EXPECT_EQ(anInter.NbSolutions(), 1); +} + +// Clear non-tangent case must keep two-line solution. +TEST_F(IntAna_IntQuadQuad_Test, CylinderCylinderNonTangent_HasTwoLines) +{ + gp_Pnt aCenter1(0.0, 0.0, 0.0); + gp_Pnt aCenter2(1.8, 0.0, 0.0); + gp_Dir aZDir(gp_Dir::D::Z); + gp_Dir anXDir(gp_Dir::D::X); + gp_Ax3 anAxis1(aCenter1, aZDir, anXDir); + gp_Ax3 anAxis2(aCenter2, aZDir, anXDir); + + gp_Cylinder aCyl1(anAxis1, 1.0); + gp_Cylinder aCyl2(anAxis2, 1.0); + + IntAna_QuadQuadGeo anInter(aCyl1, aCyl2, 1.0e-7); + + EXPECT_TRUE(anInter.IsDone()); + EXPECT_EQ(anInter.TypeInter(), IntAna_Line); + EXPECT_EQ(anInter.NbSolutions(), 2); +} + +// Sanity regression for the pre-existing skew external-tangent branch. +TEST_F(IntAna_IntQuadQuad_Test, CylinderCylinderSkewExternallyTangent_HasPoint) +{ + const double aR1 = 3.0; + const double aR2 = 2.0; + const double aDist = aR1 + aR2; + + gp_Ax3 anAxis1(gp_Pnt(0.0, 0.0, 0.0), gp_Dir::D::Z, gp_Dir::D::X); + gp_Ax3 anAxis2(gp_Pnt(0.0, aDist, 0.0), gp_Dir::D::X, gp_Dir::D::Y); + + gp_Cylinder aCyl1(anAxis1, aR1); + gp_Cylinder aCyl2(anAxis2, aR2); + + IntAna_QuadQuadGeo anInter(aCyl1, aCyl2, 1.0e-7); + + ASSERT_TRUE(anInter.IsDone()); + ASSERT_EQ(anInter.TypeInter(), IntAna_Point); + ASSERT_EQ(anInter.NbSolutions(), 1); + + const gp_Pnt aP = anInter.Point(1); + EXPECT_NEAR(gp_Lin(aCyl1.Axis()).Distance(aP), aR1, 1.0e-7); + EXPECT_NEAR(gp_Lin(aCyl2.Axis()).Distance(aP), aR2, 1.0e-7); +}