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.
This commit is contained in:
Pasukhin Dmitry
2026-04-28 12:50:12 +01:00
committed by GitHub
parent 048a2927fb
commit bea50c2fd8
6 changed files with 585 additions and 214 deletions

View File

@@ -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
)

View File

@@ -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 <gtest/gtest.h>
#include <vector>
#include <BRep_Builder.hxx>
#include <BRepAlgoAPI_Common.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
#include <BRepGProp_Face.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeWire.hxx>
#include <BRepTools.hxx>
#include <Geom_Circle.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_Plane.hxx>
#include <gp_Ax2.hxx>
#include <gp_Ax3.hxx>
#include <gp_Circ.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Pln.hxx>
#include <IntTools_Curve.hxx>
#include <IntTools_FaceFace.hxx>
#include <NCollection_Sequence.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Wire.hxx>
//==================================================================================================
// 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<IntTools_Curve>& aCurves = anInters.Lines();
EXPECT_GE(aCurves.Length() + anInters.Points().Length(), 1)
<< "Expected at least one intersection result (curve or point)";
}

View File

@@ -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
)

View File

@@ -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 <gtest/gtest.h>
#include <BRep_Builder.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_Plane.hxx>
#include <IntTools_Curve.hxx>
#include <IntTools_FaceFace.hxx>
#include <IntTools_PntOn2Faces.hxx>
#include <NCollection_Sequence.hxx>
#include <Precision.hxx>
#include <TopoDS_Face.hxx>
#include <gp_Ax3.hxx>
#include <gp_Dir.hxx>
#include <gp_Pnt.hxx>
// 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<IntTools_Curve>& aCurves = anInters.Lines();
EXPECT_GE(aCurves.Length() + anInters.Points().Length(), 1)
<< "Expected at least one intersection result (curve or point)";
}

View File

@@ -51,6 +51,7 @@
#include <gp_Cylinder.hxx>
#include <gp_Pln.hxx>
#include <gp_Sphere.hxx>
#include <MathRoot_Brent.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
@@ -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<IntPatch_WLine>& 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<IntPatch_WLine>& 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<IntPatch_WLine>& 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<IntPatch_WLine>& 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.)

View File

@@ -12,12 +12,14 @@
// commercial license or contractual agreement.
#include <IntAna_IntQuadQuad.hxx>
#include <IntAna_QuadQuadGeo.hxx>
#include <IntAna_Quadric.hxx>
#include <gp_Sphere.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Cone.hxx>
#include <gp_Pln.hxx>
#include <gp_Ax3.hxx>
#include <gp_Lin.hxx>
#include <gtest/gtest.h>
@@ -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);
}