[bos #30432] Fix intersection with edges treatment: use parameters on polygon and on triangle's edge

Added the function Tolerance2d().
Now 2d-tolerance is used instead of 3d-tolerance where it's necessary.

adm/cmake/version.cmake

@@ -19,4 +19,4 @@
 set (OCC_VERSION_MAJOR 7 )
 set (OCC_VERSION_MINOR 9 )
 set (OCC_VERSION_MAINTENANCE 0 )
-set (OCC_VERSION_DEVELOPMENT )
+set (OCC_VERSION_DEVELOPMENT "OCCT790-Salome")

src/BRep/BRep_Tool.cxx

@@ -18,6 +18,7 @@
 #include <BRep_CurveOnClosedSurface.hxx>
 #include <BRep_CurveOnSurface.hxx>
 #include <BRep_CurveRepresentation.hxx>
+#include <BRep_GCurve.hxx>
 #include <BRep_Polygon3D.hxx>
 #include <BRep_PolygonOnSurface.hxx>
 #include <BRep_PolygonOnTriangulation.hxx>
@@ -25,6 +26,7 @@
 #include <BRep_TFace.hxx>
 #include <BRep_Tool.hxx>
 #include <BRep_TVertex.hxx>
+#include <BRepAdaptor_Surface.hxx>
 #include <ElSLib.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
@@ -60,7 +62,6 @@
 #include <TopoDS_Vertex.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopTools_ShapeMapHasher.hxx>
-#include <BRep_GCurve.hxx>
 
 // modified by NIZNHY-PKV Fri Oct 17 14:13:29 2008f
 static Standard_Boolean IsPlane(const Handle(Geom_Surface)& aS);
@@ -1724,3 +1725,13 @@ Standard_Real BRep_Tool::MaxTolerance(const TopoDS_Shape&    theShape,
 
   return aTol;
 }
+
+//=======================================================================
+//function : Tolerance2d
+//purpose  : 
+//=======================================================================
+Standard_Real BRep_Tool::Tolerance2d(const TopoDS_Face& theFace, const Standard_Real theTolerance)
+{
+  BRepAdaptor_Surface aAdaptorSurface(theFace);
+  return Max(aAdaptorSurface.UResolution(theTolerance), aAdaptorSurface.VResolution(theTolerance));
+}

src/BRep/BRep_Tool.hxx

@@ -381,6 +381,13 @@ public:
   //@param theSubShape - Search subshape, only Face, Edge or Vertex are supported.
   Standard_EXPORT static Standard_Real MaxTolerance(const TopoDS_Shape&    theShape,
                                                     const TopAbs_ShapeEnum theSubShape);
+
+  //! Returns the 2d tolerance.
+  //@param theFace      - the input face to find 2d tolerance
+  //@param theTolerance - the input tolerance to calculate 2d tolerance.
+  Standard_EXPORT static Standard_Real Tolerance2d(const TopoDS_Face&  theFace,
+                                                   const Standard_Real theTolerance);
+
 };
 
 #endif // _BRep_Tool_HeaderFile

src/BRepExtrema/BRepExtrema_DistanceSS.cxx

@@ -462,9 +462,8 @@ static void PERFORM_C0(const TopoDS_Edge&         S1,
       Standard_Real Udeb, Ufin;
       BRep_Tool::Range(Eother, Udeb, Ufin);
 
-      gp_Pnt                   P1, Pt;
-      Standard_Integer         i, ii;
-      BRepClass_FaceClassifier classifier;
+      gp_Pnt           P1, Pt;
+      Standard_Integer i, ii;
       for (i = 1; i <= arrInter.Length(); i++)
       {
         const Standard_Real aParameter = arrInter(i);
@@ -796,7 +795,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Vertex&       theS1,
       Standard_Real            U, V;
       gp_Pnt                   Pt, P1 = BRep_Tool::Pnt(theS1);
       BRepClass_FaceClassifier classifier;
-      const Standard_Real      tol = BRep_Tool::Tolerance(theS2);
+      const Standard_Real      tol2d = BRep_Tool::Tolerance2d(theS2, BRep_Tool::Tolerance(theS2));
 
       for (i = 1; i <= NbExtrema; i++)
       {
@@ -808,7 +807,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Vertex&       theS1,
             // Check if the parameter does not correspond to a vertex
             Ext.Parameter(i, U, V);
             const gp_Pnt2d PUV(U, V);
-            classifier.Perform(theS2, PUV, tol);
+            classifier.Perform(theS2, PUV, tol2d);
             if (classifier.State() == TopAbs_IN)
             {
               if (myDstRef > Dstmin)
@@ -947,7 +946,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Edge&         theS1,
     if ((Dstmin < myDstRef - myEps) || (fabs(Dstmin - myDstRef) < myEps))
     {
       Standard_Real       U, V;
-      const Standard_Real tol = BRep_Tool::Tolerance(theS2);
+      const Standard_Real tol2d = BRep_Tool::Tolerance2d(theS2, BRep_Tool::Tolerance(theS2));
 
       gp_Pnt                  Pt1, Pt2;
       constexpr Standard_Real epsP = Precision::PConfusion();
@@ -966,7 +965,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Edge&         theS1,
             {
               Ext.ParameterOnFace(i, U, V);
               const gp_Pnt2d PUV(U, V);
-              classifier.Perform(theS2, PUV, tol);
+              classifier.Perform(theS2, PUV, tol2d);
               if (classifier.State() == TopAbs_IN)
               {
                 if (myDstRef > Dstmin)
@@ -999,7 +998,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Edge&         theS1,
 
     gp_Pnt              Pt;
     Standard_Real       U, V;
-    const Standard_Real tol = BRep_Tool::Tolerance(theS2);
+    const Standard_Real tol2d = BRep_Tool::Tolerance2d(theS2, BRep_Tool::Tolerance(theS2));
 
     Standard_Integer i;
     for (i = 1; i <= arrInter.Length(); i++)
@@ -1032,7 +1031,7 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Edge&         theS1,
               // Check if the parameter does not correspond to a vertex
               ExtPF.Parameter(ii, U, V);
               const gp_Pnt2d PUV(U, V);
-              classifier.Perform(theS2, PUV, tol);
+              classifier.Perform(theS2, PUV, tol2d);
               if (classifier.State() == TopAbs_IN)
               {
                 if (myDstRef > Dstmin)
@@ -1096,8 +1095,8 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Face&         theS1,
     Dstmin = sqrt(Dstmin);
     if ((Dstmin < myDstRef - myEps) || (fabs(Dstmin - myDstRef) < myEps))
     {
-      const Standard_Real tol1 = BRep_Tool::Tolerance(theS1);
-      const Standard_Real tol2 = BRep_Tool::Tolerance(theS2);
+      const Standard_Real tol2d1 = BRep_Tool::Tolerance2d(theS1, BRep_Tool::Tolerance(theS1));
+      const Standard_Real tol2d2 = BRep_Tool::Tolerance2d(theS2, BRep_Tool::Tolerance(theS2));
 
       gp_Pnt                   Pt1, Pt2;
       gp_Pnt2d                 PUV;
@@ -1115,12 +1114,12 @@ void BRepExtrema_DistanceSS::Perform(const TopoDS_Face&         theS1,
             // Check if the parameter does not correspond to a vertex
             Ext.ParameterOnFace1(i, U1, V1);
             PUV.SetCoord(U1, V1);
-            classifier.Perform(theS1, PUV, tol1);
+            classifier.Perform(theS1, PUV, tol2d1);
             if (classifier.State() == TopAbs_IN)
             {
               Ext.ParameterOnFace2(i, U2, V2);
               PUV.SetCoord(U2, V2);
-              classifier.Perform(theS2, PUV, tol2);
+              classifier.Perform(theS2, PUV, tol2d2);
               if (classifier.State() == TopAbs_IN)
               {
                 if (myDstRef > Dstmin)

src/BRepExtrema/BRepExtrema_ExtFF.cxx

@@ -80,7 +80,7 @@ void BRepExtrema_ExtFF::Perform(const TopoDS_Face& F1, const TopoDS_Face& F2)
   {
     // Exploration of points and classification
     BRepClass_FaceClassifier classifier;
-    const Standard_Real      Tol2 = BRep_Tool::Tolerance(F2);
+    const Standard_Real      Tol2d2 = BRep_Tool::Tolerance2d(F2, BRep_Tool::Tolerance(F2));
     Extrema_POnSurf          P1, P2;
 
     Standard_Integer i;
@@ -89,13 +89,14 @@ void BRepExtrema_ExtFF::Perform(const TopoDS_Face& F1, const TopoDS_Face& F2)
       myExtSS.Points(i, P1, P2);
       P1.Parameter(U1, U2);
       const gp_Pnt2d Puv1(U1, U2);
-      classifier.Perform(F1, Puv1, Tol1);
+      const Standard_Real Tol2d1 = BRep_Tool::Tolerance2d(F1, Tol1);
+      classifier.Perform(F1, Puv1, Tol2d1);
       const TopAbs_State state1 = classifier.State();
       if (state1 == TopAbs_ON || state1 == TopAbs_IN)
       {
         P2.Parameter(U1, U2);
         const gp_Pnt2d Puv2(U1, U2);
-        classifier.Perform(F2, Puv2, Tol2);
+        classifier.Perform(F2, Puv2, Tol2d2);
         const TopAbs_State state2 = classifier.State();
         if (state2 == TopAbs_ON || state2 == TopAbs_IN)
         {

src/BRepExtrema/BRepExtrema_ExtPF.cxx

@@ -79,12 +79,12 @@ void BRepExtrema_ExtPF::Perform(const TopoDS_Vertex& TheVertex, const TopoDS_Fac
   {
     BRepClass_FaceClassifier classifier;
     Standard_Real            U1, U2;
-    const Standard_Real      Tol = BRep_Tool::Tolerance(TheFace);
+    const Standard_Real      Tol2d = BRep_Tool::Tolerance2d(TheFace, BRep_Tool::Tolerance(TheFace));
     for (Standard_Integer i = 1; i <= myExtPS.NbExt(); i++)
     {
       myExtPS.Point(i).Parameter(U1, U2);
       const gp_Pnt2d Puv(U1, U2);
-      classifier.Perform(TheFace, Puv, Tol);
+      classifier.Perform(TheFace, Puv, Tol2d);
       const TopAbs_State state = classifier.State();
       if (state == TopAbs_ON || state == TopAbs_IN)
       {

src/BRepExtrema/BRepExtrema_ProximityValueTool.cxx

@@ -555,43 +555,11 @@ void BRepExtrema_ProximityValueTool::Perform(Standard_Real& theTolerance)
   if (aProximityDist1 < 0.)
     return;
 
-  // max(min) dist from the 2nd shape to t he 1st one
-  BVH_Vec3d      aP2_1, aP2_2;
-  ProxPnt_Status aPointStatus2_1 = ProxPnt_Status::ProxPnt_Status_UNKNOWN;
-  ProxPnt_Status aPointStatus2_2 = ProxPnt_Status::ProxPnt_Status_UNKNOWN;
-
-  Standard_Real aProximityDist2 = computeProximityDist(mySet2,
-                                                       myNbSamples2,
-                                                       myAddVertices2,
-                                                       myAddStatus2,
-                                                       mySet1,
-                                                       myShapeList2,
-                                                       myShapeList1,
-                                                       aP2_2,
-                                                       aP2_1,
-                                                       aPointStatus2_2,
-                                                       aPointStatus2_1);
-
-  if (aProximityDist2 < 0.)
-    return;
-
-  // min dist of the two max(min) dists
-  if (aProximityDist1 < aProximityDist2)
-  {
-    myDistance = aProximityDist1;
-    myPnt1.SetCoord(aP1_1.x(), aP1_1.y(), aP1_1.z());
-    myPnt2.SetCoord(aP1_2.x(), aP1_2.y(), aP1_2.z());
-    myPntStatus1 = aPointStatus1_1;
-    myPntStatus2 = aPointStatus1_2;
-  }
-  else
-  {
-    myDistance = aProximityDist2;
-    myPnt1.SetCoord(aP2_1.x(), aP2_1.y(), aP2_1.z());
-    myPnt2.SetCoord(aP2_2.x(), aP2_2.y(), aP2_2.z());
-    myPntStatus1 = aPointStatus2_1;
-    myPntStatus2 = aPointStatus2_2;
-  }
+  myDistance = aProximityDist1;
+  myPnt1.SetCoord(aP1_1.x(), aP1_1.y(), aP1_1.z());
+  myPnt2.SetCoord(aP1_2.x(), aP1_2.y(), aP1_2.z());
+  myPntStatus1 = aPointStatus1_1;
+  myPntStatus2 = aPointStatus1_2;
 
   myIsDone     = Standard_True;
   theTolerance = myDistance;
@@ -615,4 +583,4 @@ NCollection_CellFilter_Action BRepExtrema_VertexInspector::Inspect(const Standar
     myIsNeedAdd = Standard_False;
 
   return CellFilter_Keep;
-}
+}

src/BRepTest/BRepTest_CurveCommands.cxx

@@ -42,9 +42,12 @@
 #include <Geom2d_Circle.hxx>
 #include <Geom2dAPI_Interpolate.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
+#include <Message.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Shape.hxx>
 #include <TopoDS_Wire.hxx>
+#include <ShapeFix_Wire.hxx>
+#include <ShapeFix_Edge.hxx>
 
 #include <DBRep.hxx>
 #include <DBRep_DrawableShape.hxx>
@@ -301,6 +304,149 @@ static Standard_Integer wire(Draw_Interpretor& di, Standard_Integer n, const cha
   return 0;
 }
 
+//=======================================================================
+// strongwire
+//=======================================================================
+static Standard_Integer strongwire(Draw_Interpretor&, Standard_Integer theArgC, const char** theArgV)
+{
+  enum StrongWireMode {
+    StrongWireMode_FixTolerance = 1,
+    StrongWireMode_Approximation = 2,
+    StrongWireMode_KeepCurveType = 3
+  };
+
+  BRep_Builder aB;
+
+  TopoDS_Wire aWire;
+  aB.MakeWire(aWire);
+
+  if (theArgC < 3)
+    return 1;
+
+  // Tolerance
+  double aTolerance = Precision::Confusion();
+
+  // mode
+  StrongWireMode aMode = StrongWireMode_KeepCurveType;
+
+  // add edges
+  for (Standard_Integer anArgIter = 2; anArgIter < theArgC; anArgIter++)
+  {
+    TCollection_AsciiString aParam(theArgV[anArgIter]);
+    if (aParam == "-t")
+    {
+      anArgIter++;
+      if (anArgIter < theArgC)
+        aTolerance = Draw::Atof(theArgV[anArgIter]);
+    }
+    else if (aParam == "-m")
+    {
+      anArgIter++;
+      if (anArgIter < theArgC)
+      {
+        if (aParam == "keepType" || Draw::Atoi(theArgV[anArgIter]) == 1)
+        {
+          aMode = StrongWireMode_KeepCurveType;
+          continue;
+        }
+        else if (aParam == "approx" || Draw::Atoi(theArgV[anArgIter]) == 2)
+        {
+          aMode = StrongWireMode_Approximation;
+          continue;
+        }
+        else if (aParam == "fixTol" || Draw::Atoi(theArgV[anArgIter]) == 3)
+        {
+          aMode = StrongWireMode_FixTolerance;
+          continue;
+        }
+      }
+    }
+    else
+    {
+      TopoDS_Shape aShape = DBRep::Get(theArgV[anArgIter]);
+      if (aShape.IsNull())
+        Standard_NullObject::Raise("Shape for wire construction is null");
+      if (aShape.ShapeType() == TopAbs_EDGE || aShape.ShapeType() == TopAbs_WIRE)
+      {
+        TopExp_Explorer anExp(aShape, TopAbs_EDGE);
+        for (; anExp.More(); anExp.Next())
+          aB.Add(aWire, TopoDS::Edge(anExp.Current()));
+      }
+      else
+        Standard_TypeMismatch::Raise("Shape for wire construction is neither an edge nor a wire");
+    }
+  }
+
+  // fix edges order
+  Handle(ShapeFix_Wire) aFW = new ShapeFix_Wire;
+  aFW->Load(aWire);
+  aFW->FixReorder();
+
+  if (aFW->StatusReorder(ShapeExtend_FAIL1))
+  {
+    Message::SendFail() << "Error: Wire construction failed: several loops detected";
+    return 1;
+  }
+  else if (aFW->StatusReorder(ShapeExtend_FAIL))
+  {
+    Message::SendFail() << "Wire construction failed";
+    return 1;
+  }
+
+  bool isClosed = false;
+  Handle(ShapeAnalysis_Wire) aSaw = aFW->Analyzer();
+  if (aSaw->CheckGap3d(1)) // between last and first edges
+  {
+    Standard_Real aDist = aSaw->MinDistance3d();
+    if (aDist < aTolerance)
+      isClosed = true;
+  }
+  aFW->ClosedWireMode() = isClosed;
+  aFW->FixConnected(aTolerance);
+
+  if (aFW->StatusConnected(ShapeExtend_FAIL))
+  {
+    Message::SendFail() << "Wire construction failed: cannot build connected wire";
+    return 1;
+  }
+
+  if (aMode == StrongWireMode_KeepCurveType)
+  {
+    aFW->FixCurves();
+    if (aFW->StatusCurves(ShapeExtend_FAIL))
+    {
+      Message::SendFail() << "Wire construction failed: cannot rebuild curves through new points";
+      return 1;
+    }
+  }
+  else if (aFW->StatusConnected(ShapeExtend_DONE3))
+  {
+    if (aMode != StrongWireMode_Approximation)
+      aFW->SetPrecision(aTolerance);
+    aFW->FixGapsByRangesMode() = Standard_True;
+    if (aFW->FixGaps3d())
+    {
+      Handle(ShapeExtend_WireData) sbwd = aFW->WireData();
+      Handle(ShapeFix_Edge) aFe = new ShapeFix_Edge;
+      for (Standard_Integer anIdx = 1; anIdx <= sbwd->NbEdges(); anIdx++)
+      {
+        TopoDS_Edge aEdge = TopoDS::Edge(sbwd->Edge(anIdx));
+        aFe->FixVertexTolerance(aEdge);
+        aFe->FixSameParameter(aEdge);
+      }
+    }
+    else if (aFW->StatusGaps3d(ShapeExtend_FAIL))
+    {
+      Message::SendFail() << "Wire construction failed: cannot fix 3d gaps";
+      return 1;
+    }
+  }
+  aWire = aFW->WireAPIMake();
+
+  DBRep::Set(theArgV[1], aWire);
+  return 0;
+}
+
 //=======================================================================
 // mkedge
 //=======================================================================
@@ -2132,6 +2278,10 @@ void BRepTest::CurveCommands(Draw_Interpretor& theCommands)
 
   theCommands.Add("wire", "wire wirename [-unsorted] e1/w1 [e2/w2 ...]", __FILE__, wire, g);
 
+  theCommands.Add("strongwire",
+    "strongwire wirename e1/w1 [e2/w2 ...] [-t tol] [-m mode(keepType/1 | approx/2 | fixTol/3)]", __FILE__,
+    strongwire, g);
+
   theCommands.Add("profile", "profile, no args to get help", __FILE__, profile, g);
 
   theCommands.Add("bsplineprof", "bsplineprof, no args to get help", __FILE__, bsplineprof, g);

src/BRepTools/BRepTools.hxx

@@ -290,7 +290,7 @@ public:
           theStream,
           Standard_True,
           Standard_False,
-          TopTools_FormatVersion_CURRENT,
+          TopTools_FormatVersion_VERSION_1,
           theProgress);
   }
 
@@ -334,7 +334,7 @@ public:
                  theFile,
                  Standard_True,
                  Standard_False,
-                 TopTools_FormatVersion_CURRENT,
+                 TopTools_FormatVersion_VERSION_1,
                  theProgress);
   }
 

src/Extrema/Extrema_ExtCC.cxx

@@ -646,6 +646,15 @@ void Extrema_ExtCC::PrepareParallelResult(const Standard_Real theUt11,
                                 Precision::Confusion(),
                                 theUt11,
                                 theUt12);
+        // if (!ExtPCir.IsDone())
+        // {
+        //   // The point is in the center of the circle => there are infinite solutions
+        //   ClearSolutions();
+        //   const Standard_Real aRadius = Extrema_CurveTool::Circle(*myC[0]).Radius();
+        //   mySqDist.Append(aRadius * aRadius);
+        //   myIsParallel = Standard_True;
+        //   break;
+        // }
         if (ExtPCir.NbExt() < 1)
         {
           continue;
@@ -702,6 +711,15 @@ void Extrema_ExtCC::PrepareParallelResult(const Standard_Real theUt11,
                                 Precision::Confusion(),
                                 theUt11,
                                 theUt12);
+        // if (!ExtPCir.IsDone())
+        // {
+        //   // The point is in the center of the circle => there are infinite solutions
+        //   ClearSolutions();
+        //   const Standard_Real aRadius = Extrema_CurveTool::Circle(*myC[0]).Radius();
+        //   mySqDist.Append(aRadius * aRadius);
+        //   myIsParallel = Standard_True;
+        //   break;
+        // }
 
         Standard_Boolean isFound = !myIsParallel;
 
@@ -809,13 +827,13 @@ void Extrema_ExtCC::PrepareResults(const Extrema_ExtElC&  AlgExt,
   if (myDone)
   {
     myIsParallel = AlgExt.IsParallel();
-    if (myIsParallel)
+    NbExt = AlgExt.NbExt();
+    if (myIsParallel && NbExt == 0)
     {
       PrepareParallelResult(Ut11, Ut12, Ut21, Ut22, AlgExt.SquareDistance());
     }
     else
     {
-      NbExt = AlgExt.NbExt();
       for (i = 1; i <= NbExt; i++)
       {
         // Verification de la validite des parametres
@@ -833,11 +851,25 @@ void Extrema_ExtCC::PrepareResults(const Extrema_ExtElC&  AlgExt,
 
         if (Extrema_CurveTool::IsPeriodic(*myC[0]))
         {
-          U = ElCLib::InPeriod(U, Ut11, Ut11 + Extrema_CurveTool::Period(*myC[0]));
+          if (Abs(U - Ut11) < Precision::Confusion())
+          {
+            U = Ut11;
+          }
+          else
+          {
+            U = ElCLib::InPeriod(U, Ut11, Ut11 + Extrema_CurveTool::Period(*myC[0]));
+          }
         }
         if (Extrema_CurveTool::IsPeriodic(*myC[1]))
         {
-          U2 = ElCLib::InPeriod(U2, Ut21, Ut21 + Extrema_CurveTool::Period(*myC[1]));
+          if (Abs(U2 - Ut21) < Precision::Confusion())
+          {
+            U2 = Ut21;
+          }
+          else
+          {
+            U2 = ElCLib::InPeriod(U2, Ut21, Ut21 + Extrema_CurveTool::Period(*myC[1]));
+          }
         }
 
         if ((U >= Ut11 - RealEpsilon()) && (U <= Ut12 + RealEpsilon())

src/Extrema/Extrema_ExtElC.cxx

@@ -436,6 +436,137 @@ Standard_Boolean Extrema_ExtElC::PlanarLineCircleExtrema(const gp_Lin&  theLin,
   return Standard_True;
 }
 
+//=================================================================================================
+
+Standard_Boolean Extrema_ExtElC::PerpendicularCirclesExtrema(const gp_Circ&      C1,
+                                                             const gp_Circ&      C2,
+                                                             const gp_Pln&       aPlc1,
+                                                             const gp_Pln&       aPlc2,
+                                                             const Standard_Real aTolD,
+                                                             const Standard_Real angPlanes)
+{
+  if (angPlanes <= M_PI_4)
+  {
+    return Standard_False; // The planes are more "parallel" than "perpendicular"
+  }
+
+  Standard_Real aTolD2 = aTolD * aTolD;
+  gp_Pnt aPc1 = C1.Location();
+  gp_Pnt aPc2 = C2.Location();
+  Standard_Real aD1 = aPlc2.SquareDistance(aPc1);
+  if (aD1 < aTolD2)
+  {
+    // The center of the first circle is on the plane of the second circle,
+    // and the center of the second circle is on the plane of the first circle.
+    // => the intersection line of the two planes is along the connection of the circles' centers.
+    // Compute the number of extremes based on the distance between the centers of the circles.
+    Standard_Real aDist = aPc1.Distance(aPc2);
+    Standard_Real aMinR = Min(C1.Radius(), C2.Radius());
+    Standard_Real aMaxR = Max(C1.Radius(), C2.Radius());
+    if (aDist >= aMaxR + aMinR - aTolD)
+    {
+      // The circles are either far apart or just touching.
+      // There is one solution for the closest points on the circles.
+      myNbExt = 1;
+      gp_Vec aVec(aPc1, aPc2);
+      aVec.Normalize();
+      gp_Pnt aP1 = aPc1.Translated(aVec * C1.Radius());
+      gp_Pnt aP2 = aPc2.Translated(-aVec * C2.Radius());
+      aDist -= aMaxR + aMinR;
+      mySqDist[0] = aDist * aDist;
+      myPoint[0][0].SetValues(ElCLib::Parameter(C1, aP1), aP1);
+      myPoint[0][1].SetValues(ElCLib::Parameter(C2, aP2), aP2);
+      myIsPar = Standard_False;
+      return Standard_True;
+    }
+
+    if (Abs(aDist - aMaxR) < aTolD)
+    {
+      // The bigger circle is going through the center of the smaller circle.
+      // There are infinite solutions on the small circle.
+      myNbExt = 1;
+      gp_Vec aVec(aPc1, aPc2);
+      aVec.Normalize();
+      gp_Pnt aP1  = aPc1.Translated(aVec * C1.Radius());
+      gp_Pnt aP2  = aPc2.Translated(-aVec * C2.Radius());
+      aDist       = aMinR;
+      mySqDist[0] = aDist * aDist;
+      myPoint[0][0].SetValues(ElCLib::Parameter(C1, aP1), aP1);
+      myPoint[0][1].SetValues(ElCLib::Parameter(C2, aP2), aP2);
+      myIsPar = Standard_True; // indicator for infinite solutions
+      return Standard_True;
+    }
+
+    if (aDist <= aMaxR - aMinR + aTolD)
+    {
+      // The smaller circle is completely inside the bigger circle.
+      // There is one solution for the closest points on the circles.
+      myNbExt = 1;
+      gp_Vec aVec(aPc1, aPc2);
+      if (aVec.Magnitude() < aTolD)
+      {
+        // Both circle centers are coincident.
+        aVec = aPlc1.Axis().Direction().Crossed(aPlc2.Axis().Direction());
+        myNbExt++;
+      }
+      aVec.Normalize();
+      if (C1.Radius() < C2.Radius())
+      {
+        aVec.Reverse();
+      }
+      gp_Pnt        aP1 = aPc1.Translated(aVec * C1.Radius());
+      gp_Pnt        aP2 = aPc2.Translated(aVec * C2.Radius());
+      Standard_Real aD  = aMaxR - aDist - aMinR;
+      mySqDist[0]       = aD * aD;
+      myPoint[0][0].SetValues(ElCLib::Parameter(C1, aP1), aP1);
+      myPoint[0][1].SetValues(ElCLib::Parameter(C2, aP2), aP2);
+      if (myNbExt > 1)
+      {
+        gp_Pnt aP3 = aPc1.Translated(-aVec * C1.Radius());
+        gp_Pnt aP4 = aPc2.Translated(-aVec * C2.Radius());
+        mySqDist[1] = mySqDist[0];
+        myPoint[1][0].SetValues(ElCLib::Parameter(C1, aP3), aP3);
+        myPoint[1][1].SetValues(ElCLib::Parameter(C2, aP4), aP4);
+      }
+      myIsPar = Standard_False;
+      return Standard_True;
+    }
+
+    if (aDist < aMaxR + aMinR - aTolD)
+    {
+      // The circles are intersecting.
+      // There is one solution for the closest points on the circles.
+      myNbExt = 1;
+      gp_Vec aVec(aPc2, aPc1);
+      if (aVec.Magnitude() < aMaxR)
+      {
+        // The center of the small circle is inside the big circle.
+        aVec.Reverse();
+      }
+      aVec.Normalize();
+      if (C1.Radius() < C2.Radius())
+      {
+        aVec.Reverse();
+      }
+      gp_Pnt        aP1 = aPc1.Translated(aVec * C1.Radius());
+      gp_Pnt        aP2 = aPc2.Translated(aVec * C2.Radius());
+      Standard_Real aD  = (aMaxR + aMinR - aDist);
+      mySqDist[0]       = aD * aD;
+      myPoint[0][0].SetValues(ElCLib::Parameter(C1, aP1), aP1);
+      myPoint[0][1].SetValues(ElCLib::Parameter(C2, aP2), aP2);
+      myIsPar = Standard_False;
+      return Standard_True;
+    }
+
+    // std::cout << "====> Extrema_ExtElC: unsupported case for circles in different planes." << std::endl;
+    // std::cout << "  Distance between centers: " << aDist << std::endl;
+    // std::cout << "  Angle between planes ...: " << angPlanes * 180.0 / M_PI << " deg" << std::endl;
+    // std::cout << "  Radius Circle 1 ........: " << C1.Radius() << std::endl;
+    // std::cout << "  Radius Circle 2 ........: " << C2.Radius() << std::endl;
+  }
+  return Standard_False;
+}
+
 //=======================================================================
 // function : Extrema_ExtElC
 // purpose  :
@@ -973,6 +1104,22 @@ Extrema_ExtElC::Extrema_ExtElC(const gp_Circ& C1, const gp_Circ& C2)
   bIsSamePlane = aDc1.IsParallel(aDc2, aTolA) && aD2 < aTolD2;
   if (!bIsSamePlane)
   {
+    // Handle the specific case:
+    //  * where both planes are "almost perpendicular" (=> angle between planes > 45 deg)
+    //  * where the center of each circle is on the plane of the other circle
+    //    (=> the circle centers are on the intersection line of the two planes)
+    const Standard_Real angPlanes = aDc1.Angle(aDc2);
+    if (aD2 < aTolD2 && angPlanes > M_PI_4)
+    {
+      // The center of the second circle is on the plane of the first circle and
+      // both planes are "almost perpendicular".
+      gp_Pln aPlc2(aPc2, aDc2);
+      if (PerpendicularCirclesExtrema(C1, C2, aPlc1, aPlc2, aTolD, angPlanes))
+      {
+        myDone = Standard_True;
+        return;
+      }
+    }
     return;
   }
 

src/Extrema/Extrema_ExtElC.hxx

@@ -28,6 +28,7 @@ class gp_Circ;
 class gp_Elips;
 class gp_Hypr;
 class gp_Parab;
+class gp_Pln;
 
 //! It calculates all the distance between two elementary
 //! curves.
@@ -85,7 +86,13 @@ public:
 protected:
   //! Computes extrema in case when considered line and circle are in one plane
   Standard_EXPORT Standard_Boolean PlanarLineCircleExtrema(const gp_Lin& C1, const gp_Circ& C2);
-
+  //! Computes extrema in case when two circles are in almost perpendicular planes
+  //! and their centers are on the intersection line of these planes
+  //! (angle between planes > 45 degrees)
+  Standard_EXPORT Standard_Boolean PerpendicularCirclesExtrema(const gp_Circ& C1, const gp_Circ& C2,
+                                                               const gp_Pln& aPlc1, const gp_Pln& aPlc2,
+                                                               const Standard_Real aTolD,
+                                                               const Standard_Real angPlanes);
 private:
   Standard_Boolean myDone;
   Standard_Boolean myIsPar;

src/Intf/Intf_InterferencePolygonPolyhedron.gxx

@@ -1010,7 +1010,7 @@ void Intf_InterferencePolygonPolyhedron::Intersect
 	      if(IsInSegment(VecPol,gp_Vec(BegO,PO),NVecPol,ParamOnO,Tolerance)) {
 		if(IsInSegment(VecTri,gp_Vec(BegT,PT),NVecTri,ParamOnT,Tolerance)) {
 		  //-- cout<<" * "<<endl;
-		  gp_XYZ spLieu=BegT.XYZ()+((EndT.XYZ()-BegT.XYZ())*param);
+		  gp_XYZ spLieu=BegO.XYZ()+((EndO.XYZ()-BegO.XYZ())*ParamOnO);
 		  Standard_Integer tmin,tmax;
 		  if(pTri_i>pTri_ip1pc3) { 
 		    tmin=pTri_ip1pc3; tmax=pTri_i; 
@@ -1022,7 +1022,7 @@ void Intf_InterferencePolygonPolyhedron::Intersect
 				       typOnG, 0, iLin, ParamOnO, 
 				       Intf_EDGE, 
 				       tmin, 
-				       tmax, 0.,
+				       tmax, ParamOnT,
 				       1.);
 		  mySPoins.Append(SP);
 		}

src/ShapeFix/ShapeFix_Wire.cxx

@@ -53,13 +53,19 @@
 #include <BRepBuilderAPI_MakeEdge.hxx>
 #include <BRepBuilderAPI_MakeVertex.hxx>
 #include <BRepTools.hxx>
+#include <GC_MakeArcOfCircle.hxx>
+#include <GC_MakeArcOfEllipse.hxx>
+#include <GC_MakeLine.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2d_Line.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
 #include <Geom2dConvert.hxx>
-#include <Geom_BSplineCurve.hxx>
 #include <Geom_Curve.hxx>
+#include <Geom_Circle.hxx>
+#include <Geom_Ellipse.hxx>
+#include <Geom_BezierCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
 #include <Geom_OffsetCurve.hxx>
 #include <Geom_Plane.hxx>
 #include <Geom_SphericalSurface.hxx>
@@ -69,7 +75,10 @@
 #include <GeomAdaptor_Curve.hxx>
 #include <GeomAPI.hxx>
 #include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
+#include <gp_Elips.hxx>
+#include <gp_Pln.hxx>
 #include <IntRes2d_IntersectionPoint.hxx>
 #include <IntRes2d_SequenceOfIntersectionPoint.hxx>
 #include <Message_Msg.hxx>
@@ -216,6 +225,7 @@ void ShapeFix_Wire::ClearStatuses()
   myStatusReorder          = emptyStatus;
   myStatusSmall            = emptyStatus;
   myStatusConnected        = emptyStatus;
+  myStatusCurves           = emptyStatus;
   myStatusEdgeCurves       = emptyStatus;
   myStatusDegenerated      = emptyStatus;
   myStatusSelfIntersection = emptyStatus;
@@ -497,6 +507,21 @@ Standard_Boolean ShapeFix_Wire::FixConnected(const Standard_Real prec)
 
 //=================================================================================================
 
+Standard_Boolean ShapeFix_Wire::FixCurves()
+{
+  myStatusCurves = ShapeExtend::EncodeStatus(ShapeExtend_OK);
+  if (!IsLoaded()) return Standard_False;
+
+  for (Standard_Integer anIdx = NbEdges(); anIdx > 0; anIdx--) {
+    FixCurves(anIdx);
+    myStatusCurves |= myLastFixStatus;
+  }
+
+  return StatusCurves(ShapeExtend_DONE);
+}
+
+//=================================================================================================
+
 Standard_Boolean ShapeFix_Wire::FixEdgeCurves()
 {
   myStatusEdgeCurves = ShapeExtend::EncodeStatus(ShapeExtend_OK);
@@ -1365,6 +1390,193 @@ Standard_Boolean ShapeFix_Wire::FixConnected(const Standard_Integer num, const S
 
 //=================================================================================================
 
+Standard_Boolean ShapeFix_Wire::FixCurves(const Standard_Integer theIdx)
+{
+  // assume fix curves step should be after fix vertices
+  myLastFixStatus = ShapeExtend::EncodeStatus(ShapeExtend_OK);
+  if (!IsLoaded() || NbEdges() <= 0) 
+    return Standard_False;
+
+  // action: replacing curves in edges
+  Handle(ShapeExtend_WireData) anSbwd = WireData();
+  Standard_Integer anIdx = (theIdx > 0 ? theIdx : anSbwd->NbEdges());
+  TopoDS_Edge anEdge = anSbwd->Edge(anIdx);
+  Standard_Real aPrec = BRep_Tool::Tolerance(anEdge);
+  ShapeAnalysis_Edge aSae;
+  ShapeAnalysis_Wire aSaw;
+  Handle(Geom_Curve) aCurve3d;
+  Standard_Real aCurBounds[3];
+  Standard_Boolean IsReversed = Standard_False;
+  aSae.Curve3d(anEdge, aCurve3d, aCurBounds[0], aCurBounds[2], IsReversed);
+  aCurBounds[1] = (aCurBounds[0] + aCurBounds[2]) / 2;
+  gp_Pnt anEnds[3];
+  anEnds[0] = BRep_Tool::Pnt(aSae.FirstVertex(anEdge));
+  aCurve3d->D0(aCurBounds[1], anEnds[1]);
+  anEnds[2] = BRep_Tool::Pnt(aSae.LastVertex(anEdge));
+
+  gp_Pnt aGeomEnds[2];
+  aCurve3d->D0(aCurBounds[0], aGeomEnds[0]);
+  aCurve3d->D0(aCurBounds[2], aGeomEnds[1]);
+
+  Standard_Real aGap0 = Min(anEnds[0].Distance(aGeomEnds[0]), anEnds[0].Distance(aGeomEnds[1]));
+  Standard_Real aGap2 = Min(anEnds[2].Distance(aGeomEnds[0]), anEnds[2].Distance(aGeomEnds[1]));
+  if (Max (aGap0, aGap2) < Precision::Confusion()) // nothing to do
+    return true;
+
+  if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Circle))) {
+    Standard_Real anOldR = Handle(Geom_Circle)::DownCast(aCurve3d)->Circ().Radius();
+    gp_Vec anArcNorm = gp_Vec(anEnds[2], anEnds[0]) / 2;
+    gp_Pnt aCenter(anEnds[0].XYZ() - anArcNorm.XYZ());
+    int aSign = anEnds[1].Distance(aCenter) > anOldR ? 1 : -1; // for arc length > PI
+    Standard_Real aD = anOldR*anOldR - anArcNorm.SquareMagnitude();
+    aD = abs(aD) < Precision::Confusion() ? 0. : aD;
+    Standard_Real anArcR = anOldR + aSign * sqrt(aD);
+    gp_Ax2 aNormal(aCenter, anArcNorm);
+    Handle(Geom_Circle) anArc = new Geom_Circle(aNormal, anArcR);
+    GeomAPI_ProjectPointOnCurve projector(anEnds[1], anArc);
+    anEnds[1] = projector.NearestPoint();
+    GC_MakeArcOfCircle arc(anEnds[0], anEnds[1], anEnds[2]);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(arc.Value()).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Ellipse))) {
+    /// aaajfa: BEGIN - to provide elliptic edge FORWARD orientation
+    gp_Pnt tmpPnt = anEnds[0];
+    anEnds[0] = anEnds[2];
+    anEnds[2] = tmpPnt;
+    /// aaajfa: END - to provide elliptic edge FORWARD orientation
+
+    Handle(Geom_Ellipse) anOld = Handle(Geom_Ellipse)::DownCast(aCurve3d);
+    Handle(Geom_Plane) aPln = new Geom_Plane(anEnds[0], gp_Vec(anEnds[2], anEnds[0]).Crossed(gp_Vec(anEnds[1], anEnds[0])));
+    GeomAPI_ProjectPointOnSurf aProjector(anOld->Elips().Location(), aPln);
+    gp_Pnt anOrigin = aProjector.NearestPoint();
+    aProjector.Init(anOld->Elips().Location().XYZ() + anOld->Elips().XAxis().Direction().XYZ(), aPln);
+    gp_Ax2 anAx(anOrigin, aPln->Axis().Direction(), aProjector.NearestPoint().XYZ() - anOrigin.XYZ());
+
+    // compute angle
+    Standard_Real aRec = DBL_MAX;
+    Standard_Real anAngle = 0.;
+    Standard_Integer aSplNum = 10;
+    for (Standard_Integer anIdxI = -aSplNum; anIdxI < aSplNum; ++anIdxI)
+    {
+      Handle(Geom_Ellipse) anEll = new Geom_Ellipse(anAx, anOld->MajorRadius(), anOld->MinorRadius());
+      Standard_Real anAnglei = aPrec*anIdxI / anEll->MajorRadius() / aSplNum;
+      anEll->Rotate(anAx.Axis(), anAnglei);
+      GeomAPI_ProjectPointOnCurve aProjector1(anEnds[0], anEll);
+      Standard_Real aDist = 0.;
+      for (Standard_Integer anIdxJ = 0; anIdxJ < 2; ++anIdxJ)
+      {
+        aProjector1.Perform(anEnds[anIdxJ *2]);
+        aDist += std::fmin (aProjector1.Distance(1), aProjector1.Distance(2));
+      }
+      if (aDist < aRec)
+      {
+        aRec = aDist;
+        anAngle = anAnglei;
+      }
+    }
+    gp_Elips aTemp(anAx, anOld->MajorRadius(), anOld->MinorRadius());
+    aTemp.Rotate(anAx.Axis(), anAngle);
+
+    // compute shift
+    gp_Vec aX = aTemp.XAxis().Direction();
+    gp_Vec aY = aTemp.YAxis().Direction();
+    gp_Pnt2d aP1((anEnds[0].XYZ() - anOrigin.XYZ()).Dot(aX.XYZ()), (anEnds[0].XYZ() - anOrigin.XYZ()).Dot(aY.XYZ()));
+    gp_Pnt2d aP2((anEnds[2].XYZ() - anOrigin.XYZ()).Dot(aX.XYZ()), (anEnds[2].XYZ() - anOrigin.XYZ()).Dot(aY.XYZ()));
+
+    // x = ky + p   linear equation
+    // where (x, y) shift point, 
+    // k, p constant coefficients
+    Standard_Real k = 1, p = 0;
+    Standard_Real R = anOld->MajorRadius();
+    Standard_Real r = anOld->MinorRadius();
+    k = (R / r) * (R / r) *
+      (aP1.Y() - aP2.Y()) / (aP2.X() - aP1.X());
+    p = -(1. / 2) * (R / r) * (R / r) *
+      (aP1.Y()*aP1.Y() - aP2.Y()*aP2.Y()) / (aP2.X() - aP1.X()) + aP1.X() / 2 + aP2.X() / 2;
+    // ax^2 + bx + c = 0  square equation 
+    // a, b, c constant coefficients
+    Standard_Real a = 0., b = 0., c = 0.;
+    a = R*R + k*k*r*r;
+    b = 2 * (k*p*r*r - k*aP1.X()*r*r - aP1.Y()*R*R);
+    c = aP1.X()*aP1.X()*r*r +
+      aP1.Y()*aP1.Y()*R*R -
+      r*r*R*R +
+      p*p*r*r - 2 * aP1.X()*p*r*r;
+    Standard_Real y1 = (-b - sqrt(b*b - 4 * a*c)) / 2 / a;
+    Standard_Real y2 = (-b + sqrt(b*b - 4 * a*c)) / 2 / a;
+    Standard_Real x1 = k*y1 + p;
+    Standard_Real x2 = k*y2 + p;
+
+    gp_Pnt anOri = anOld->Location();
+    if (x1*x1 + y1*y1 < x2*x2 + y2*y2) 
+      anOri = anOri.XYZ() + aX.XYZ()*x1 + aY.XYZ()*y1;
+    else
+      anOri = anOri.XYZ() + aX.XYZ()*x2 + aY.XYZ()*y2;
+    aTemp.SetLocation(anOri);
+
+    GC_MakeArcOfEllipse anArc(aTemp, anEnds[2], anEnds[0], true);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(anArc.Value()).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_Line)))
+  {
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(anEnds[0], anEnds[2]).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) 
+  {
+    Handle(Geom_BSplineCurve) anOld = Handle(Geom_BSplineCurve)::DownCast(aCurve3d);
+
+    if (anOld->Pole(1).Distance(aGeomEnds[0]) > Precision::Confusion() ||
+        anOld->Pole(anOld->NbPoles()).Distance(aGeomEnds[1]) > Precision::Confusion()) {
+      // FAIL1 means we cannot fix Bezier or B-Spline curve
+      // because its ends do not correspond to first and last poles
+      // (i.e. it is a piece of entire curve)
+      myLastFixStatus |= ShapeExtend::EncodeStatus ( ShapeExtend_FAIL1 );
+      return false;
+    }
+
+    Handle(Geom_Geometry) aNewG = anOld->Copy();
+    Handle(Geom_BSplineCurve) aNewC = Handle(Geom_BSplineCurve)::DownCast(aNewG);
+    int p = anEnds[0].Distance(aGeomEnds[0]) < anEnds[1].Distance(aGeomEnds[0]) ? 0 : 2;
+    aNewC->SetPole(1, anEnds[p]);
+    aNewC->SetPole(anOld->NbPoles(), anEnds[2-p]);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(aNewC).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+  else if (aCurve3d->IsKind(STANDARD_TYPE(Geom_BezierCurve))) 
+  {
+    Handle(Geom_BezierCurve) anOld = Handle(Geom_BezierCurve)::DownCast(aCurve3d);
+
+    if (anOld->Pole(1).Distance(aGeomEnds[0]) > Precision::Confusion() ||
+        anOld->Pole(anOld->NbPoles()).Distance(aGeomEnds[1]) > Precision::Confusion()) {
+      // FAIL1 means we cannot fix Bezier or B-Spline curve
+      // because its ends do not correspond to first and last poles
+      // (i.e. it is a piece of entire curve)
+      myLastFixStatus |= ShapeExtend::EncodeStatus ( ShapeExtend_FAIL1 );
+      return false;
+    }
+
+    Handle(Geom_Geometry) aNewG = anOld->Copy();
+    Handle(Geom_BezierCurve) aNewC = Handle(Geom_BezierCurve)::DownCast(aNewG);
+    int p = anEnds[0].Distance(aGeomEnds[0]) < anEnds[1].Distance(aGeomEnds[0]) ? 0 : 2;
+    aNewC->SetPole(1, anEnds[p]);
+    aNewC->SetPole(anOld->NbPoles(), anEnds[2-p]);
+    TopoDS_Edge aNewEdge = BRepBuilderAPI_MakeEdge(aNewC).Edge();
+    anSbwd->Set(aNewEdge, theIdx);
+    return true;
+  }
+
+  return true;
+}
+
+//=================================================================================================
+
 Standard_Boolean ShapeFix_Wire::FixSeam(const Standard_Integer num)
 {
   myLastFixStatus = ShapeExtend::EncodeStatus(ShapeExtend_OK);

src/ShapeFix/ShapeFix_Wire.hxx

@@ -280,6 +280,9 @@ public:
   //! If <prec> is -1 then MaxTolerance() is taken.
   Standard_EXPORT Standard_Boolean FixConnected(const Standard_Real prec = -1.0);
 
+  //! Applies FixCurves(num) to all edges in the wire
+  Standard_EXPORT Standard_Boolean FixCurves();
+  
   //! Groups the fixes dealing with 3d and pcurves of the edges.
   //! The order of the fixes and the default behaviour are:
   //! ShapeFix_Edge::FixReversed2d
@@ -352,6 +355,9 @@ public:
   Standard_EXPORT Standard_Boolean FixConnected(const Standard_Integer num,
                                                 const Standard_Real    prec);
 
+  //! Fixes curves
+  Standard_EXPORT Standard_Boolean FixCurves(const Standard_Integer num);
+  
   //! Fixes a seam edge
   //! A Seam edge has two pcurves, one for forward. one for reversed
   //! The forward pcurve must be set as first
@@ -420,6 +426,8 @@ public:
 
   Standard_Boolean StatusConnected(const ShapeExtend_Status status) const;
 
+  Standard_Boolean StatusCurves(const ShapeExtend_Status status) const;
+
   Standard_Boolean StatusEdgeCurves(const ShapeExtend_Status status) const;
 
   Standard_Boolean StatusDegenerated(const ShapeExtend_Status status) const;
@@ -500,6 +508,7 @@ protected:
   Standard_Integer           myStatusReorder;
   Standard_Integer           myStatusSmall;
   Standard_Integer           myStatusConnected;
+  Standard_Integer           myStatusCurves;
   Standard_Integer           myStatusEdgeCurves;
   Standard_Integer           myStatusDegenerated;
   Standard_Integer           myStatusClosed;

src/ShapeFix/ShapeFix_Wire.lxx

@@ -447,6 +447,16 @@ inline Standard_Boolean ShapeFix_Wire::StatusConnected(const ShapeExtend_Status
   return ShapeExtend::DecodeStatus(myStatusConnected, status);
 }
 
+//=======================================================================
+// function : StatusCurves
+// purpose  :
+//=======================================================================
+
+inline  Standard_Boolean ShapeFix_Wire::StatusCurves(const ShapeExtend_Status status) const
+{
+  return ShapeExtend::DecodeStatus(myStatusCurves, status); 
+}
+
 //=======================================================================
 // function : StatusEdgeCurves
 // purpose  :

src/StdSelect/StdSelect_BRepSelectionTool.cxx

@@ -66,45 +66,10 @@
 
 #define BVH_PRIMITIVE_LIMIT 800000
 
-namespace
-{
-//=======================================================================
-// function : getCylinderCircles
-// purpose  : Extracts up to two circular edges from a hollow cylinder face
-//=======================================================================
-std::array<gp_Circ, 2> getCylinderCircles(const TopoDS_Face& theHollowCylinder,
-                                          Standard_Size&     theNumCircles)
-{
-  std::array<gp_Circ, 2> aCircles; // Array to store up to two circles
-  theNumCircles             = 0;   // Initialize the number of circles found
-  Standard_Integer aLinesNb = 0;   // Counter for the number of edges processed
-
-  for (TopExp_Explorer anEdgeExp(theHollowCylinder, TopAbs_EDGE); anEdgeExp.More();
-       anEdgeExp.Next())
-  {
-    const TopoDS_Edge& anEdge = TopoDS::Edge(anEdgeExp.Current());
-    BRepAdaptor_Curve  anAdaptor(anEdge);
-    aLinesNb++;
-
-    if (anAdaptor.GetType() == GeomAbs_Circle && BRep_Tool::IsClosed(anEdge))
-    {
-      theNumCircles++;
-      aCircles[theNumCircles - 1] = anAdaptor.Circle();
-    }
-    else if (anAdaptor.GetType() != GeomAbs_Line || aLinesNb > 4)
-    {
-      theNumCircles = 0;
-      return std::array<gp_Circ, 2>();
-    }
-    if (theNumCircles == 2)
-    {
-      break;
-    }
-  }
-
-  return aCircles;
-}
-} // namespace
+// The following flag can be used to enable optimization of cone/cylinder selection
+// Unfortunately, this optimization is not stable and may lead to incorrect selection
+// in some cases. It is disabled by default.
+//#define OPTIMIZE_CONE_CYLINDER_SELECTION
 
 //==================================================
 // function: PreBuildBVH
@@ -642,6 +607,48 @@ void StdSelect_BRepSelectionTool::GetEdgeSensitive(const TopoDS_Shape&
   }
 }
 
+#ifdef OPTIMIZE_CONE_CYLINDER_SELECTION
+namespace
+{
+//=======================================================================
+// function : getCylinderCircles
+// purpose  : Extracts up to two circular edges from a hollow cylinder face
+//=======================================================================
+std::array<gp_Circ, 2> getCylinderCircles(const TopoDS_Face& theHollowCylinder,
+                                          Standard_Size&     theNumCircles)
+{
+  std::array<gp_Circ, 2> aCircles; // Array to store up to two circles
+  theNumCircles             = 0;   // Initialize the number of circles found
+  Standard_Integer aLinesNb = 0;   // Counter for the number of edges processed
+
+  for (TopExp_Explorer anEdgeExp(theHollowCylinder, TopAbs_EDGE); anEdgeExp.More();
+       anEdgeExp.Next())
+  {
+    const TopoDS_Edge& anEdge = TopoDS::Edge(anEdgeExp.Current());
+    BRepAdaptor_Curve  anAdaptor(anEdge);
+    aLinesNb++;
+
+    if (anAdaptor.GetType() == GeomAbs_Circle && BRep_Tool::IsClosed(anEdge))
+    {
+      theNumCircles++;
+      aCircles[theNumCircles - 1] = anAdaptor.Circle();
+    }
+    else if (anAdaptor.GetType() != GeomAbs_Line || aLinesNb > 4)
+    {
+      theNumCircles = 0;
+      return std::array<gp_Circ, 2>();
+    }
+    if (theNumCircles == 2)
+    {
+      break;
+    }
+  }
+
+  return aCircles;
+}
+} // namespace
+#endif /*OPTIMIZE_CONE_CYLINDER_SELECTION*/
+
 //=================================================================================================
 
 Standard_Boolean StdSelect_BRepSelectionTool::GetSensitiveForFace(
@@ -688,6 +695,7 @@ Standard_Boolean StdSelect_BRepSelectionTool::GetSensitiveForFace(
         return Standard_True;
       }
     }
+#ifdef OPTIMIZE_CONE_CYLINDER_SELECTION
     else if (Handle(Geom_ConicalSurface) aGeomCone = Handle(Geom_ConicalSurface)::DownCast(aSurf))
     {
       Standard_Size                aNumCircles;
@@ -748,6 +756,7 @@ Standard_Boolean StdSelect_BRepSelectionTool::GetSensitiveForFace(
         return Standard_True;
       }
     }
+#endif /*OPTIMIZE_CONE_CYLINDER_SELECTION*/
     else if (Handle(Geom_Plane) aGeomPlane = Handle(Geom_Plane)::DownCast(aSurf))
     {
       TopTools_IndexedMapOfShape aSubfacesMap;

tests/bugs/modalg_8/bug33394

@@ -0,0 +1,14 @@
+puts "========================"
+puts "0033394: Modeling Algorithms - Wrong usage of BRepClass_FaceClassifier with 3d tolerance"
+puts "========================"
+puts ""
+
+pload MODELING
+restore [locate_data_file bug33394.brep] fb
+vertex vert 1016.3915670000133 17180 -5689.3758362036951
+distmini aa fb vert
+set tolerance [checkmaxtol aa]
+
+if { [dval aa_val] > ${tolerance} } {
+  puts "Error: The distance should be 0"
+}

tests/heal/grids.list

@@ -25,4 +25,5 @@
 025 update_tolerance_locked
 026 checkshape
 027 split_number
-028 split_two_numbers
+028 split_two_numbers
+029 wire_fix_curves

tests/heal/wire_fix_curves/strongwire_circle

@@ -0,0 +1,16 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+circle c1 0 0 0 15
+circle c2 0 0 0 15
+trim c1 c1 0 pi
+trim c2 c2 pi 2*pi
+translate c1 0 0 0.00000005
+translate c2 0 0 -0.00000005
+mkedge c1 c1
+mkedge c2 c2
+strongwire result c1 c2
+checkshape result
+whatis result

tests/heal/wire_fix_curves/strongwire_ellipse

@@ -0,0 +1,16 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+ellipse e1 0 0 0 25 15
+ellipse e2 0 0 0 25 15
+trim e1 e1 0 pi
+trim e2 e2 pi 2*pi
+translate e1 0 0 0.00000005
+translate e2 0 0 -0.00000005
+mkedge e1 e1
+mkedge e2 e2
+strongwire result e1 e2
+checkshape result
+whatis result

tests/heal/wire_fix_curves/strongwire_halfcircle1

@@ -0,0 +1,24 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+point px 100   0   0
+point py   0 100   0
+point pz   0   0 100
+
+vertex vx px
+vertex vy py
+
+edge exy vx vy
+
+gcarc arc cir py pz px
+
+# Here we have an error during automatic grid execution with test and testgrid:
+# Tcl Exception: EMPTY
+# But in line-by-line entering or script sourcing mode this line doesn't produce any error
+
+mkedge earc arc
+
+strongwire result exy earc
+whatis result

tests/heal/wire_fix_curves/strongwire_halfcircle2

@@ -0,0 +1,25 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+point p1   0   0  0
+point p2 100   0  0
+point p4 100 110  0
+
+gcarc arc cir p1 p2 p4
+
+# Here we have an error during automatic grid execution with test and testgrid:
+# Tcl Exception: EMPTY
+# But in line-by-line entering or script sourcing mode this line doesn't produce any error
+
+mkedge arc1 arc
+mkedge arc2 arc
+
+tmirror arc2 0 0 0 -110 100 0 -copy
+ttranslate arc2 0 0 0.003 -copy
+
+strongwire result arc1 arc2 -t 0.01 -m keepType
+
+checkshape result
+whatis result

tests/heal/wire_fix_curves/strongwire_nurb

@@ -0,0 +1,14 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+bsplinecurve c1 3 2  -1.0 4   1.0 4  0 0 0 1  1 4 0 1  2 4 0 1  3 0 0 1
+bsplinecurve c2 3 2  -1.0 4   1.0 4  0 0 0 1  1 -4 0 1  2 -4 0 1  3 0 0 1
+translate c1 0 0 0.00000005
+translate c2 0 0 -0.00000005
+mkedge c1 c1
+mkedge c2 c2
+strongwire result c1 c2
+checkshape result
+whatis result

tests/heal/wire_fix_curves/strongwire_triangle

@@ -0,0 +1,17 @@
+#############################################################
+## Fix wire algo
+## moves a curve of invalid wire to make it valid
+#############################################################
+
+vertex v1   0   0 0
+vertex v2 100   0 0
+vertex v3  50 100 0
+
+edge e1 v1 v2
+edge e2 v2 v3
+edge e3 v3 v1
+
+strongwire result e1 e2 e3
+
+checkshape result
+whatis result