From 4d60da8c58839e9a64741be513cdebb2acd9eab0 Mon Sep 17 00:00:00 2001
From: ifv <ifv@opencascade.com>
Date: Tue, 6 Dec 2022 10:12:22 +0300
Subject: [PATCH] 0033244: Modeling Algorithms - Surface-surface intersection
produces the double curves
IntAna_QuadQuadGeo.cxx - estimation of angular tolerance is added for case cone-cone
tests/lowalgos/intss/bug33244 - new test case added
---
src/IntAna/IntAna_QuadQuadGeo.cxx | 123 ++++++++++++++++++++++++++----
tests/lowalgos/intss/bug33244 | 48 ++++++++++++
2 files changed, 157 insertions(+), 14 deletions(-)
create mode 100644 tests/lowalgos/intss/bug33244
diff --git a/src/IntAna/IntAna_QuadQuadGeo.cxx b/src/IntAna/IntAna_QuadQuadGeo.cxx
index 3ab50bcd85..c3dcf2c767 100644
--- a/src/IntAna/IntAna_QuadQuadGeo.cxx
+++ b/src/IntAna/IntAna_QuadQuadGeo.cxx
@@ -50,11 +50,21 @@
#include <Standard_DomainError.hxx>
#include <Standard_OutOfRange.hxx>
#include <StdFail_NotDone.hxx>
+#include <gce_MakePln.hxx>
+#include <ProjLib.hxx>
+#include <IntAna2d_AnaIntersection.hxx>
+#include <IntAna2d_IntPoint.hxx>
+
+#ifdef DEBUGLINES
+#include <Geom2d_Line.hxx>
+#endif
static
gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
static
void RefineDir(gp_Dir& aDir);
+static
+ Standard_Real EstimDist(const gp_Cone& theCon1, const gp_Cone& theCon2);
//=======================================================================
//class : AxeOperator
@@ -249,6 +259,72 @@ gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
}
}
+
+//=======================================================================
+//function : EstimDist
+//purpose : returns a minimal distance from apex to any solution
+//=======================================================================
+Standard_Real EstimDist(const gp_Cone& theCon1, const gp_Cone& theCon2)
+{
+ //It is supposed that axes of cones are coplanar and
+ //distance between them > Precision::Confusion()
+ gp_Pnt aPA1 = theCon1.Apex(), aPA2 = theCon2.Apex();
+
+ gp_Pnt aP3 = aPA1.Translated(theCon1.Position().Direction());
+
+ gce_MakePln aMkPln(aPA1, aPA2, aP3);
+ if(!aMkPln.IsDone())
+ return Precision::Infinite();
+
+ const gp_Pln& aPln = aMkPln.Value();
+
+ gp_Lin anAx1(aPA1, theCon1.Position().Direction());
+ gp_Lin2d anAx12d = ProjLib::Project(aPln, anAx1);
+ gp_Lin2d Lines1[2];
+ Standard_Real anAng1 = theCon1.SemiAngle();
+ Lines1[0] = anAx12d.Rotated(anAx12d.Location(), anAng1);
+ Lines1[1] = anAx12d.Rotated(anAx12d.Location(), -anAng1);
+ //
+ gp_Lin anAx2(aPA2, theCon2.Position().Direction());
+ gp_Lin2d anAx22d = ProjLib::Project(aPln, anAx2);
+ gp_Lin2d Lines2[2];
+ Standard_Real anAng2 = theCon2.SemiAngle();
+ Lines2[0] = anAx22d.Rotated(anAx22d.Location(), anAng2);
+ Lines2[1] = anAx22d.Rotated(anAx22d.Location(), -anAng2);
+
+#ifdef DEBUGLINES
+ Handle(Geom2d_Line) L10 = new Geom2d_Line(Lines1[0]);
+ Handle(Geom2d_Line) L11 = new Geom2d_Line(Lines1[1]);
+ Handle(Geom2d_Line) L20 = new Geom2d_Line(Lines2[0]);
+ Handle(Geom2d_Line) L21 = new Geom2d_Line(Lines2[1]);
+#endif
+
+ Standard_Real aMinDist[2] = { Precision::Infinite(), Precision::Infinite() };
+ Standard_Integer i, j, k;
+ IntAna2d_AnaIntersection anInter;
+ for (i = 0; i < 2; ++i)
+ {
+ for (j = 0; j < 2; ++j)
+ {
+ anInter.Perform(Lines1[i], Lines2[j]);
+ if (anInter.IsDone())
+ {
+ Standard_Integer aNbPoints = anInter.NbPoints();
+ for (k = 1; k <= aNbPoints; ++k)
+ {
+ const IntAna2d_IntPoint& anIntP = anInter.Point(k);
+ Standard_Real aPar1 = Abs(anIntP.ParamOnFirst());
+ aMinDist[0] = Min(aPar1, aMinDist[0]);
+ Standard_Real aPar2 = Abs(anIntP.ParamOnSecond());
+ aMinDist[1] = Min(aPar2, aMinDist[1]);
+ }
+ }
+ }
+ }
+
+ Standard_Real aDist = Max(aMinDist[0], aMinDist[1]);
+ return aDist;
+}
//=======================================================================
//function : IntAna_QuadQuadGeo
//purpose : Empty constructor
@@ -1356,30 +1432,49 @@ IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
//=======================================================================
void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
const gp_Cone& Con2,
- const Standard_Real Tol)
+ const Standard_Real Tol)
{
- done=Standard_True;
+ done = Standard_True;
//
Standard_Real tg1, tg2, aDA1A2, aTol2;
gp_Pnt aPApex1, aPApex2;
Standard_Real TOL_APEX_CONF = 1.e-10;
-
- //
- tg1=Tan(Con1.SemiAngle());
- tg2=Tan(Con2.SemiAngle());
- if((tg1 * tg2) < 0.) {
+ //
+ tg1 = Tan(Con1.SemiAngle());
+ tg2 = Tan(Con2.SemiAngle());
+
+ if ((tg1 * tg2) < 0.) {
tg2 = -tg2;
}
//
- aTol2=Tol*Tol;
- aPApex1=Con1.Apex();
- aPApex2=Con2.Apex();
- aDA1A2=aPApex1.SquareDistance(aPApex2);
+ aTol2 = Tol*Tol;
+ aPApex1 = Con1.Apex();
+ aPApex2 = Con2.Apex();
+ aDA1A2 = aPApex1.SquareDistance(aPApex2);
//
- AxeOperator A1A2(Con1.Axis(),Con2.Axis());
+ AxeOperator A1A2(Con1.Axis(), Con2.Axis());
//
+ Standard_Real aTolAng = myEPSILON_ANGLE_CONE;
+ if ((Abs(tg1 - tg2) < Tol) && (A1A2.Parallel()))
+ {
+ Standard_Real DistA1A2 = A1A2.Distance();
+ if (DistA1A2 > 100. * Tol)
+ {
+ Standard_Real aMinSolDist = EstimDist(Con1, Con2);
+ if (aMinSolDist < Epsilon(1.))
+ {
+ aTolAng = Tol;
+ }
+ else
+ {
+ aTolAng = Max(myEPSILON_ANGLE_CONE, Tol / aMinSolDist);
+ aTolAng = Min(aTolAng, Tol);
+ }
+ }
+ }
+
// 1
if(A1A2.Same()) {
//-- two circles
@@ -1427,8 +1522,8 @@ IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
}
} //-- fin A1A2.Same
// 2
- else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
- //-- voir AnVer12mai98
+ else if((Abs(tg1-tg2) < aTolAng ) && (A1A2.Parallel())) {
+
Standard_Real DistA1A2=A1A2.Distance();
gp_Dir DA1=Con1.Position().Direction();
gp_Vec O1O2(Con1.Apex(),Con2.Apex());
diff --git a/tests/lowalgos/intss/bug33244 b/tests/lowalgos/intss/bug33244
new file mode 100644
index 0000000000..8a571926a0
--- /dev/null
+++ b/tests/lowalgos/intss/bug33244
@@ -0,0 +1,48 @@
+
+puts "============"
+puts "0033244: Modeling Algorithms - Surface-surface intersection produces the double curves"
+puts "============"
+puts ""
+
+
+set PI180 0.017453292519943295
+
+set x -1.11630646267172
+set y -4.54487349779333
+set z 13.2493435203532
+set dx -1.05794851588922e-07
+set dy -1.39278337794573e-08
+set dz 0.999999999999994
+set ux 0.999999999999994
+set uy 5.91645678915759e-31
+set uz 1.05794851588922e-07
+set semi-angle [expr 0.785398163360967 / ${PI180}]
+set radius 0.560000000061149
+
+cone s1 ${x} ${y} ${z} ${dx} ${dy} ${dz} ${semi-angle} ${radius}
+
+set x -2.08647872350287e-07
+set y -5.78732475509323
+set z 13.2493436211
+set dx -1.05794850062242e-07
+set dy -1.39278350756825e-08
+set dz 0.999999999999994
+set ux 0.999999999999995
+set uy 0
+set uz 1.05794850062242e-07
+set semi-angle [expr 0.785398163396248 / ${PI180}]
+set radius 0.785398163396248
+
+cone s2 ${x} ${y} ${z} ${dx} ${dy} ${dz} ${semi-angle} ${radius}
+
+intersect ii s1 s2
+if { ![isdraw ii_1] || ![isdraw ii_2] || [isdraw ii_3] } {
+ puts "ERROR. Intersection is wrong"
+}
+
+intersect jj s2 s1
+if { ![isdraw jj_1] || ![isdraw jj_2] || [isdraw jj_3] } {
+ puts "ERROR. Intersection is wrong"
+}
+
+