0029807: [Regression to 7.0.0] Impossible to cut cone from prism

The algorithm has been improved for the cases when the intersection line goes through the cone apex.

<!break>

1. All special points are put to the ALine forcefully (if they are true intersection point). Currently this step has not been implemented yet.

2. Now the tolerance of IntPatch_Point (put into ALine) is computed in order to cover the distance between it and the correspond ALine.

3. Test cases have been created.

4. Procedure of trimming IntAna_Curve has been improved.

5. Criterion when the discriminant of IntAna_Curve can be considered to be equal to 0 has been improved.

6. Methods IntAna_Curve::FindParameter(...) (and IntPatch_ALine::FindParameter(...)) currently returns list of all parameters corresponding the given point (IntAna_Curve can be self-interfered curve). Before the fix, this method always returned only one (randomly chosen) parameter.

7. Interfaces of the following methods have been changed: IntAna_Curve::FindParameter(...), IntPatch_ALine::FindParameter(...), IntPatch_ALine::ChangeVertex(...), IntPatch_SpecialPoints::AddPointOnUorVIso(...), IntPatch_SpecialPoints::AddSingularPole(...), IntPatch_WLineTool::ExtendTwoWLines().

8. Following methods have been added: IntAna_Quadric::SpecialPoints(...), IntPatch_ALineToWLine::GetSectionRadius(...), IntPatch_SpecialPoints::ProcessSphere(...), IntPatch_SpecialPoints::ProcessCone(...), IntPatch_SpecialPoints::GetTangentToIntLineForCone(...).

------------------
1) tests/boolean/volumemaker/C5
   tests/boolean/volumemaker/C6
   tests/boolean/volumemaker/E7

They are real IMPROVEMENTS. In the FIX (in compare with MASTER), section result between pairs of faces f2&f6 (C5), f3&f7 (C6) and f1&f5 (E7) is closed. Separated test cases have been created in order to focus on the problem with section. Bug #28503 has been fixed.

Correction in test cases.

src/BRepTest/BRepTest_CheckCommands.cxx

@@ -298,9 +298,25 @@ static Standard_Integer checksection(Draw_Interpretor& di,
 				     Standard_Integer narg, const char** a)
 {
   if (narg < 2) {
+    di << a[0] << " shape [-r <ref_val>]\n";
     return 1;
   }
+
+  Standard_Integer aCompareValue = -1;
   TopoDS_Shape S = DBRep::Get(a[1]);
+
+  for (Standard_Integer anAI = 2; anAI < narg; anAI++)
+  {
+    if (!strcmp(a[anAI], "-r"))
+    {
+      aCompareValue = Draw::Atoi(a[++anAI]);
+    }
+    else
+    {
+      di << "Error: Wrong option" << a[anAI] << "\n";
+    }
+  }
+
   TopTools_MapOfShape theVertices;
   TopExp_Explorer exp;
   for (exp.Init(S, TopAbs_VERTEX); exp.More(); exp.Next()) {
@@ -309,6 +325,20 @@ static Standard_Integer checksection(Draw_Interpretor& di,
   }
   //cout << " nb alone Vertices : " << theVertices.Extent() << endl;
   di << " nb alone Vertices : " << theVertices.Extent() << "\n";
+
+  if (aCompareValue >= 0)
+  {
+    if (theVertices.Extent() == aCompareValue)
+    {
+      di << "Section is OK\n";
+    }
+    else
+    {
+      di << "Error: "<< aCompareValue << " vertices are expected but " <<
+                        theVertices.Extent() << " are found.\n";
+    }
+  }
+
   char Name[32];
   Standard_Integer ipp=0;
   TopTools_MapIteratorOfMapOfShape itvx;
@@ -1698,7 +1728,8 @@ void BRepTest::CheckCommands(Draw_Interpretor& theCommands)
 //  Modified by skv - Tue Apr 27 13:35:39 2004 End
 
   theCommands.Add("checksection", 
-		  "checks the closure of a section : checksection name",
+		  "checks the closure of a section : checksection name [-r <RefVal>]\n"
+                  "\"-r\" - allowed number of allone vertices.",
 		  __FILE__,
 		  checksection,
 		  g);

src/Geom/Geom_ConicalSurface.hxx

@@ -40,7 +40,7 @@ class Geom_ConicalSurface;
 DEFINE_STANDARD_HANDLE(Geom_ConicalSurface, Geom_ElementarySurface)
 
 //! Describes a cone.
-//! A cone is defined by the half-angle at its apex, and
+//! A cone is defined by the half-angle (can be negative) at its apex, and
 //! is positioned in space by a coordinate system (a
 //! gp_Ax3 object) and a reference radius as follows:
 //! - The "main Axis" of the coordinate system is the
@@ -79,7 +79,8 @@ public:
   
 
   //! A3 defines the local coordinate system of the conical surface.
-  //! Ang is the conical surface semi-angle ]0, PI/2[.
+  //! Ang is the conical surface semi-angle. Its absolute value is in range
+  //! ]0, PI/2[.
   //! Radius is the radius of the circle Viso in the placement plane
   //! of the conical surface defined with "XAxis" and "YAxis".
   //! The "ZDirection" of A3 defines the direction of the surface's
@@ -90,8 +91,8 @@ public:
   //! such that the normal Vector (N = D1U ^ D1V) is oriented towards
   //! the "outside region" of the surface.
   //!
-  //! Raised if Radius < 0.0 or Ang < Resolution from gp or
+  //! Raised if Radius < 0.0 or Abs(Ang) < Resolution from gp or
-  //! Ang >= PI/2 - Resolution
+  //! Abs(Ang) >= PI/2 - Resolution
   Standard_EXPORT Geom_ConicalSurface(const gp_Ax3& A3, const Standard_Real Ang, const Standard_Real Radius);
   
 
@@ -112,9 +113,11 @@ public:
   
 
   //! Changes the semi angle of the conical surface.
-  //!
+  //! Semi-angle can be negative. Its absolute value
-  //! Raised if Ang < Resolution or Ang >= PI/2 - Resolution
+  //! Abs(Ang) is in range ]0,PI/2[.
-  Standard_EXPORT void SetSemiAngle (const Standard_Real Ang);
+  //! Raises ConstructionError if Abs(Ang) < Resolution from gp or
+  //! Abs(Ang) >= PI/2 - Resolution
+  Standard_EXPORT void SetSemiAngle(const Standard_Real Ang);
   
 
   //! returns a non transient cone with the same geometric properties
@@ -206,7 +209,8 @@ public:
   Standard_EXPORT Standard_Real RefRadius() const;
   
 
-  //! returns the semi-angle of the conical surface ]0.0, PI/2[.
+  //! Returns the semi-angle at the apex of this cone.
+  //! Attention! Semi-angle can be negative.
   Standard_EXPORT Standard_Real SemiAngle() const;
   
   //! returns True.

src/IntAna/IntAna_Curve.cxx

@@ -37,6 +37,8 @@
 //--       pas etre mene a bien.  
 //----------------------------------------------------------------------
 
+#include <algorithm>
+
 #include <ElSLib.hxx>
 #include <gp_Cone.hxx>
 #include <gp_Cylinder.hxx>
@@ -52,7 +54,7 @@
 //function : IntAna_Curve
 //purpose  : 
 //=======================================================================
-  IntAna_Curve::IntAna_Curve()
+IntAna_Curve::IntAna_Curve()
 {
   typequadric=GeomAbs_OtherSurface;
   firstbounded=Standard_False;
@@ -62,22 +64,22 @@
 //function : SetConeQuadValues
 //purpose  : Description de l intersection Cone Quadrique
 //=======================================================================
-  void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
+void IntAna_Curve::SetConeQuadValues(const gp_Cone& Cone,
-				       const Standard_Real Qxx,
+                                     const Standard_Real Qxx,
-				       const Standard_Real Qyy,
+                                     const Standard_Real Qyy,
-				       const Standard_Real Qzz,
+                                     const Standard_Real Qzz,
-				       const Standard_Real Qxy,
+                                     const Standard_Real Qxy,
-				       const Standard_Real Qxz,
+                                     const Standard_Real Qxz,
-				       const Standard_Real Qyz,
+                                     const Standard_Real Qyz,
-				       const Standard_Real Qx,
+                                     const Standard_Real Qx,
-				       const Standard_Real Qy,
+                                     const Standard_Real Qy,
-				       const Standard_Real Qz,
+                                     const Standard_Real Qz,
-				       const Standard_Real Q1,
+                                     const Standard_Real Q1,
-				       const Standard_Real TOL,
+                                     const Standard_Real TOL,
-				       const Standard_Real DomInf,
+                                     const Standard_Real DomInf,
-				       const Standard_Real DomSup,
+                                     const Standard_Real DomSup,
-				       const Standard_Boolean twocurves,
+                                     const Standard_Boolean twocurves,
-				       const Standard_Boolean takezpositive) 
+                                     const Standard_Boolean takezpositive)
 {
   
   Ax3        = Cone.Position();
@@ -112,36 +114,40 @@
   Z2Cos      = (UnSurTgAngle+UnSurTgAngle)*Qxz;
   Z2CosCos   = Qxx;
   Z2SinSin   = Qyy;
-  Z2CosSin   = Qxy+Qxy;
+  Z2CosSin   = Qxy;
 
   Tolerance  = TOL;
-  DomainInf  = DomInf;
+  DomainInf = DomInf;
-  DomainSup  = DomSup;
+  DomainSup = DomSup;
   
   RestrictedInf = RestrictedSup = Standard_True;   //-- Le Domaine est Borne
   firstbounded = lastbounded = Standard_False;
+
+  myFirstParameter = DomainInf;
+  myLastParameter = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
+                                  DomainSup;
 }
 
 //=======================================================================
 //function : SetCylinderQuadValues
 //purpose  : Description de l intersection Cylindre Quadrique
 //=======================================================================
-  void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
+void IntAna_Curve::SetCylinderQuadValues(const gp_Cylinder& Cyl,
-					   const Standard_Real Qxx,
+                                         const Standard_Real Qxx,
-					   const Standard_Real Qyy,
+                                         const Standard_Real Qyy,
-					   const Standard_Real Qzz,
+                                         const Standard_Real Qzz,
-					   const Standard_Real Qxy,
+                                         const Standard_Real Qxy,
-					   const Standard_Real Qxz,
+                                         const Standard_Real Qxz,
-					   const Standard_Real Qyz,
+                                         const Standard_Real Qyz,
-					   const Standard_Real Qx,
+                                         const Standard_Real Qx,
-					   const Standard_Real Qy,
+                                         const Standard_Real Qy,
-					   const Standard_Real Qz,
+                                         const Standard_Real Qz,
-					   const Standard_Real Q1,
+                                         const Standard_Real Q1,
-					   const Standard_Real TOL,
+                                         const Standard_Real TOL,
-					   const Standard_Real DomInf,
+                                         const Standard_Real DomInf,
-					   const Standard_Real DomSup,
+                                         const Standard_Real DomSup,
-					   const Standard_Boolean twocurves,
+                                         const Standard_Boolean twocurves,
-					   const Standard_Boolean takezpositive) 
+                                         const Standard_Boolean takezpositive)
 {
   
   Ax3        = Cyl.Position();
@@ -157,7 +163,7 @@
   Z0Cos      = RCylmul2*Qx;
   Z0CosCos   = Qxx*RCyl*RCyl;
   Z0SinSin   = Qyy*RCyl*RCyl;
-  Z0CosSin   = RCylmul2*RCyl*Qxy;
+  Z0CosSin   = RCyl*RCyl*Qxy;
 
   Z1Cte      = Qz+Qz;
   Z1Sin      = RCylmul2*Qyz;
@@ -174,18 +180,22 @@
   Z2CosSin   = 0.0;
 
   Tolerance  = TOL;
-  DomainInf  = DomInf;
+  DomainInf = DomInf;
-  DomainSup  = DomSup;
+  DomainSup = DomSup;
   
   RestrictedInf = RestrictedSup = Standard_True;
   firstbounded = lastbounded = Standard_False;
+
+  myFirstParameter = DomainInf;
+  myLastParameter  = (TwoCurves) ? DomainSup + DomainSup - DomainInf :
+                                   DomainSup;
 }
 
 //=======================================================================
 //function : IsOpen
 //purpose  : 
 //=======================================================================
-  Standard_Boolean IntAna_Curve::IsOpen() const 
+Standard_Boolean IntAna_Curve::IsOpen() const
 {
   return(RestrictedInf && RestrictedSup);
 }
@@ -194,17 +204,16 @@
 //function : Domain
 //purpose  : 
 //=======================================================================
-  void IntAna_Curve::Domain(Standard_Real& DInf,
+void IntAna_Curve::Domain(Standard_Real& theFirst,
-			    Standard_Real& DSup) const
+                          Standard_Real& theLast) const
 {
-  if(RestrictedInf && RestrictedSup) {
+  if (RestrictedInf && RestrictedSup)
-    DInf=DomainInf;
+  {
-    DSup=DomainSup;
+    theFirst = myFirstParameter;
-    if(TwoCurves) {
+    theLast = myLastParameter;
-      DSup+=DSup-DInf;
-    }
   }
-  else {
+  else
+  {
     throw Standard_DomainError("IntAna_Curve::Domain");
   }
 }
@@ -212,7 +221,7 @@
 //function : IsConstant
 //purpose  : 
 //=======================================================================
-  Standard_Boolean IntAna_Curve::IsConstant() const 
+Standard_Boolean IntAna_Curve::IsConstant() const
 {
   //-- ???  Pas facile de decider a la seule vue des Param.
   return(Standard_False);
@@ -222,7 +231,7 @@
 //function : IsFirstOpen
 //purpose  : 
 //=======================================================================
-  Standard_Boolean IntAna_Curve::IsFirstOpen() const 
+Standard_Boolean IntAna_Curve::IsFirstOpen() const
 {
   return(firstbounded);
 }
@@ -231,7 +240,7 @@
 //function : IsLastOpen
 //purpose  : 
 //=======================================================================
-  Standard_Boolean IntAna_Curve::IsLastOpen() const 
+Standard_Boolean IntAna_Curve::IsLastOpen() const
 {
   return(lastbounded);
 }
@@ -239,7 +248,7 @@
 //function : SetIsFirstOpen
 //purpose  : 
 //=======================================================================
-  void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag) 
+void IntAna_Curve::SetIsFirstOpen(const Standard_Boolean Flag)
 {
   firstbounded = Flag;
 }
@@ -248,7 +257,7 @@
 //function : SetIsLastOpen
 //purpose  : 
 //=======================================================================
-  void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag) 
+void IntAna_Curve::SetIsLastOpen(const Standard_Boolean Flag)
 {
   lastbounded = Flag;
 }
@@ -257,29 +266,40 @@
 //function : InternalUVValue
 //purpose  : 
 //=======================================================================
-  void IntAna_Curve::InternalUVValue(const Standard_Real theta,
+void IntAna_Curve::InternalUVValue(const Standard_Real theta,
-				     Standard_Real& Param1,
+                                   Standard_Real& Param1,
-				     Standard_Real& Param2,
+                                   Standard_Real& Param2,
-				     Standard_Real& A,
+                                   Standard_Real& A,
-				     Standard_Real& B,
+                                   Standard_Real& B,
-				     Standard_Real& C,
+                                   Standard_Real& C,
-				     Standard_Real& cost,
+                                   Standard_Real& cost,
-				     Standard_Real& sint,
+                                   Standard_Real& sint,
-				     Standard_Real& SigneSqrtDis) const 
+                                   Standard_Real& SigneSqrtDis) const
 {
   const Standard_Real aRelTolp = 1.0+Epsilon(1.0), aRelTolm = 1.0-Epsilon(1.0);
+  
+  // Infinitesimal step of increasing curve parameter. See comment below.
+  const Standard_Real aDT = 100.0*Epsilon(DomainSup + DomainSup - DomainInf);
+
   Standard_Real Theta=theta;
   Standard_Boolean SecondSolution=Standard_False; 
 
-  if((Theta<DomainInf*aRelTolm)  ||  
+  if ((Theta<DomainInf*aRelTolm) ||
-     ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||  
+      ((Theta>DomainSup*aRelTolp) && (!TwoCurves)) ||
-     (Theta>(DomainSup+DomainSup-DomainInf)*aRelTolp)) {
+      (Theta>(DomainSup + DomainSup - DomainInf)*aRelTolp))
+  {
     SigneSqrtDis = 0.;
     throw Standard_DomainError("IntAna_Curve::Domain");
   }
   
-  if(Theta>DomainSup) { 
+  if (Abs(Theta - DomainSup) < aDT)
-    Theta=DomainSup+DomainSup-Theta;
+  {
+    // Point of Null-discriminant.
+    Theta = DomainSup;
+  }
+  else if (Theta>DomainSup)
+  {
+    Theta = DomainSup + DomainSup - Theta;
     SecondSolution=Standard_True; 
   }
 
@@ -291,53 +311,56 @@
   //
   cost = Cos(Theta);
   sint = Sin(Theta);
-  Standard_Real costsint = cost*sint;
+  const Standard_Real aSin2t = Sin(Theta + Theta);
+  const Standard_Real aCos2t = Cos(Theta + Theta);
   
   A=Z2Cte+sint*(Z2Sin+sint*Z2SinSin)+cost*(Z2Cos+cost*Z2CosCos)
-    +Z2CosSin*costsint;
+    + Z2CosSin*aSin2t;
+  
+  const Standard_Real aDA = cost*Z2Sin - sint*Z2Cos + 
+                            aSin2t*(Z2SinSin - Z2CosCos) + 
+                            aCos2t*(Z2CosSin * Z2CosSin);
 
   B=Z1Cte+sint*(Z1Sin+sint*Z1SinSin)+cost*(Z1Cos+cost*Z1CosCos)
-    +Z1CosSin*costsint;
+    + Z1CosSin*aSin2t;
+
+  const Standard_Real aDB = Z1Sin*cost - Z1Cos*sint +
+                            aSin2t*(Z1SinSin - Z1CosCos) +
+                            aCos2t*(Z1CosSin + Z1CosSin);
   
   C=Z0Cte+sint*(Z0Sin+sint*Z0SinSin)+cost*(Z0Cos+cost*Z0CosCos)
-    +Z0CosSin*costsint;
+    + Z0CosSin*aSin2t;
   
+  const Standard_Real aDC = Z0Sin*cost - Z0Cos*sint +
+                            aSin2t*(Z0SinSin - Z0CosCos) + 
+                            aCos2t*(Z0CosSin + Z0CosSin);
 
-  const Standard_Real aDiscriminant = Max(B*B-4.0*A*C, 0.0);
+  Standard_Real aDiscriminant = B*B-4.0*A*C;
 
-  if(Abs(A)<=Precision::PConfusion()) {   
+  // We consider that infinitesimal dt = aDT.
-    //-- cout<<" IntAna_Curve:: Internal UV Value : A="<<A<<"  -> Abs(A)="<<Abs(A)<<endl;
+  // Error of discriminant computation is equal to
-    if(Abs(B)<=Precision::PConfusion()) { 
+  // (d(Disc)/dt)*dt, where 1st derivative d(Disc)/dt = 2*B*aDB - 4*(A*aDC + C*aDA).
-      //-- cout<<" Probleme :  Pas de solutions "<<endl;
+
-      Param2=0.0;
+  const Standard_Real aTolD = 2.0*aDT*Abs(B*aDB - 2.0*(A*aDC + C*aDA));
+  
+  if (aDiscriminant < aTolD)
+    aDiscriminant = 0.0;
+
+  if (Abs(A) <= Precision::PConfusion())
+  {
+    if (Abs(B) <= Precision::PConfusion())
+    {
+      Param2 = 0.0;
     }
-    else {
+    else
-      //modified by NIZNHY-PKV Fri Dec  2 16:02:46 2005f
+    {
-      Param2 = -C/B;
+      Param2 = -C / B;
-      /*
-      if(!SecondSolution) {   
-	//-- Cas Param2 = (-B+Sqrt(Discriminant))/(A+A);
-	//--            = (-B+Sqrt(B**2 - Eps)) / 2A
-	//--            = -C / B
-	Param2 = -C/B;
-      } 
-      else {
-	//-- Cas Param2 = (-B-Sqrt(Discriminant))/(A+A);
-	//--            = (-B-Sqrt(B**2 - Eps)) / 2A
-	if(A) { 
-	  Param2 = -B/A; 
-	}
-	else { 
-	  Param2 = -B*10000000.0;
-	}
-      }
-      */
-      //modified by NIZNHY-PKV Fri Dec  2 16:02:54 2005t
     }
   }
-  else {
+  else
-    SigneSqrtDis = (SecondSolution)? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
+  {
-    Param2=(-B+SigneSqrtDis)/(A+A);
+    SigneSqrtDis = (SecondSolution) ? Sqrt(aDiscriminant) : -Sqrt(aDiscriminant);
+    Param2 = (-B + SigneSqrtDis) / (A + A);
   }
 }
 
@@ -345,7 +368,7 @@
 //function : Value
 //purpose  : 
 //=======================================================================
-  gp_Pnt IntAna_Curve::Value(const Standard_Real theta)  
+gp_Pnt IntAna_Curve::Value(const Standard_Real theta)
 {
   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
   //
@@ -361,9 +384,9 @@
 //function : D1u
 //purpose  : 
 //=======================================================================
-  Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
+Standard_Boolean IntAna_Curve::D1u(const Standard_Real theta,
-				     gp_Pnt& Pt,
+                                   gp_Pnt& Pt,
-				     gp_Vec& Vec) 
+                                   gp_Vec& Vec)
 {
   //-- Pour detecter le cas ou le calcul est impossible 
   Standard_Real A, B, C, U, V, sint, cost, SigneSqrtDis;
@@ -374,14 +397,15 @@
   InternalUVValue(theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
   //
   Pt = Value(theta);
-  if(Abs(A)<0.0000001 || Abs(SigneSqrtDis)<0.0000000001) return(Standard_False);
+  if (Abs(A)<1.0e-7 || Abs(SigneSqrtDis)<1.0e-10) return(Standard_False);
   
   
   //-- Approximation de la derivee (mieux que le calcul mathematique!)
-  Standard_Real dtheta = (DomainSup-DomainInf)*0.000001;
+  Standard_Real dtheta = (DomainSup - DomainInf)*1.0e-6;
   Standard_Real theta2 = theta+dtheta;
-  if((theta2<DomainInf)  ||  ((theta2>DomainSup) && (!TwoCurves))
+  if ((theta2<DomainInf) || ((theta2>DomainSup) && (!TwoCurves))
-     ||  (theta2>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
+      || (theta2>(DomainSup + DomainSup - DomainInf + 1.0e-14)))
+  {
     dtheta = -dtheta;
     theta2 = theta+dtheta;
   }
@@ -395,147 +419,93 @@
 }  
 //=======================================================================
 //function : FindParameter
-//purpose  : Para est en sortie le parametre sur la courbe 
+//purpose  : Projects P to the ALine. Returns the list of parameters as a results
+//            of projection.
+//           Sometimes aline can be self-intersected line (see bug #29807 where 
+//            ALine goes through the cone apex).
 //=======================================================================
-  Standard_Boolean IntAna_Curve::FindParameter (const gp_Pnt& P,
+void IntAna_Curve::FindParameter(const gp_Pnt& theP,
-						Standard_Real& Para) const
+                                 TColStd_ListOfReal& theParams) const
 {
-  Standard_Real theta,z, aTolPrecision=0.0001;
+  const Standard_Real aPIpPI = M_PI + M_PI,
-  Standard_Real PIpPI = M_PI + M_PI;
+                      anEpsAng = 1.e-8,
+                      aSqTolPrecision=1.0e-8;
+  Standard_Real aTheta = 0.0;
   //
-  switch (typequadric) {
+  switch (typequadric)
-
+  {
-  case GeomAbs_Cylinder:
+    case GeomAbs_Cylinder:
     {
-      ElSLib::CylinderParameters(Ax3,RCyl,P,theta,z);
+      Standard_Real aZ;
+      ElSLib::CylinderParameters(Ax3, RCyl, theP, aTheta, aZ);
     }
     break;
 
-  case GeomAbs_Cone :
+    case GeomAbs_Cone:
     {
-      ElSLib::ConeParameters(Ax3,RCyl,Angle,P,theta,z);
+      Standard_Real aZ;
+      ElSLib::ConeParameters(Ax3, RCyl, Angle, theP, aTheta, aZ);
     }
     break;
 
-  default:
+    default:
-    return Standard_False;
+      return;
-    break;
   }
   //
-  Standard_Real epsAng = 1.e-8;
+  if (!firstbounded && (DomainInf > aTheta) && ((DomainInf - aTheta) <= anEpsAng))
-  Standard_Real tmin = DomainInf;
+  {
-  Standard_Real tmax = DomainSup;
+    aTheta = DomainInf;
-  Standard_Real U,V,A,B,C,sint,cost,SigneSqrtDis;
-  Standard_Real z1,z2;
-
-  A=0.0;  B=0.0;   C=0.0;
-  U=0.0;  V=0.0;  
-  sint=0.0;  cost=0.0;
-  SigneSqrtDis=0.0;
-  //U=V=A=B=C=sint=cost=SigneSqrtDis=0.0;
-  //
-  if (!firstbounded && tmin > theta && (tmin-theta) <= epsAng) {
-    theta = tmin;
   }
-  else if (!lastbounded && theta > tmax && (theta-tmax) <= epsAng) {
+  else if (!lastbounded && (aTheta > DomainSup) && ((aTheta - DomainSup) <= anEpsAng))
-    theta = tmax;
+  {
+    aTheta = DomainSup;
   }
   //
-  if (theta < tmin ) {
+  if (aTheta < DomainInf)
-    theta = theta + PIpPI;
+  {
+    aTheta = aTheta + aPIpPI;
   }
-  else if (theta > tmax) {
+  else if (aTheta > DomainSup)
-    theta = theta - PIpPI;
+  {
-  }
+    aTheta = aTheta - aPIpPI;
-  if (theta < tmin || theta > tmax) {
-    if(theta>tmax) { 
-      InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis); 
-      gp_Pnt PMax(InternalValue(U,V)); 
-      if(PMax.Distance(P) < aTolPrecision) { 
-	Para = tmax; 
-	return(Standard_True); 
-      }
-    }
-    if(theta<tmin) { 
-      InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis); 
-      gp_Pnt PMin(InternalValue(U,V));  
-      if(PMin.Distance(P) < aTolPrecision) { 
-	Para = tmin; 
-	return(Standard_True); 
-      }
-    }
-    //-- lbr le 14 Fev 96 : On teste malgre tout si le point n est pas le 
-    //-- point de debut ou de fin 
-    //-- cout<<"False 1 "<<endl;
-    //    theta = tmin;                               le 25 Nov 96
   }
 
-  if (TwoCurves) {
+  const Standard_Integer aMaxPar = 5;
-    if(theta > tmax) 
+  Standard_Real aParams[aMaxPar] = {DomainInf, DomainSup, aTheta,
-      theta =  tmax;
+                                    (TwoCurves)? DomainSup + DomainSup - aTheta : RealLast(),
-    if(theta < tmin) 
+                                    (TwoCurves) ? DomainSup + DomainSup - DomainInf : RealLast()};
-      theta = tmin;
-    InternalUVValue(theta,U,z1,A,B,C,cost,sint,SigneSqrtDis); 
-    A = B = C = sint = cost = SigneSqrtDis = 0.0;
-    InternalUVValue(tmax+tmax - theta,U,z2,A,B,C,cost,sint,SigneSqrtDis); 
 
-    if (Abs(z-z1) <= Abs(z-z2)) {
+  std::sort(aParams, aParams + aMaxPar - 1);
-      Para = theta;
+
-    }
+  for (Standard_Integer i = 0; i < aMaxPar; i++)
-    else {
+  {
-      Para = tmax+tmax - theta;
+    if (aParams[i] > myLastParameter)
+      break;
+    
+    if (aParams[i] < myFirstParameter)
+      continue;
+
+    if (i && (aParams[i] - aParams[i - 1]) < Precision::PConfusion())
+      continue;
+
+    Standard_Real U = 0.0, V= 0.0, 
+                  A = 0.0, B = 0.0, C = 0.0,
+                  sint = 0.0, cost = 0.0, SigneSqrtDis = 0.0;
+    InternalUVValue(aParams[i], U, V, A, B, C, 
+                    cost, sint, SigneSqrtDis);
+    const gp_Pnt aP(InternalValue(U, V));
+    if (aP.SquareDistance(theP) < aSqTolPrecision)
+    {
+      theParams.Append(aParams[i]);
     }
   }
-  else {
-    Para = theta;
-  }
-
-  if((Para<DomainInf)  ||  ((Para>DomainSup) && (!TwoCurves))
-     ||  (Para>(DomainSup+DomainSup-DomainInf+0.00000000000001))) {
-    return(Standard_False);
-  }
-  
-  InternalUVValue(Para,U,V,A,B,C,cost,sint,SigneSqrtDis);
-  gp_Pnt PPara = InternalValue(U,V);
-  Standard_Real Dist = PPara.Distance(P);
-  if(Dist > aTolPrecision) { 
-    //-- Il y a eu un probleme  
-    //-- On teste si le point est un point double 
-    InternalUVValue(tmin,U,V,A,B,C,cost,sint,SigneSqrtDis);
-    PPara = InternalValue(U,V);
-    Dist = PPara.Distance(P);
-    if(Dist <= aTolPrecision) {
-      Para = tmin; 
-      return(Standard_True);
-    }
-
-    InternalUVValue(tmax,U,V,A,B,C,cost,sint,SigneSqrtDis);
-    PPara = InternalValue(U,V);
-    Dist = PPara.Distance(P);
-    if(Dist <= aTolPrecision) {
-      Para = tmax; 
-      return(Standard_True);
-    }
-    if (TwoCurves) {
-      Standard_Real Theta = DomainSup+DomainSup-DomainInf;
-      InternalUVValue(Theta,U,V,A,B,C,cost,sint,SigneSqrtDis);
-      PPara = InternalValue(U,V);
-      Dist = PPara.Distance(P);
-      if(Dist <= aTolPrecision) {
-	Para = Theta; 
-	return(Standard_True);
-      } 
-    }
-    return(Standard_False);
-  }
-  return(Standard_True);
 }
 //=======================================================================
 //function : InternalValue
 //purpose  : 
 //=======================================================================
-  gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
+gp_Pnt IntAna_Curve::InternalValue(const Standard_Real U,
-				     const Standard_Real _V) const 
+                                   const Standard_Real _V) const
 {
   //-- cout<<" ["<<U<<","<<V<<"]";
   Standard_Real V = _V;
@@ -569,13 +539,14 @@
 //function : SetDomain
 //purpose  : 
 //=======================================================================
-  void IntAna_Curve::SetDomain(const Standard_Real DInf,
+void IntAna_Curve::SetDomain(const Standard_Real theFirst,
-			       const Standard_Real DSup) 
+                             const Standard_Real theLast)
 {
-  if(DInf>=DSup) {
+  if (theLast <= theFirst)
+  {
     throw Standard_DomainError("IntAna_Curve::Domain");
   }
   //
-  DomainInf=DInf;
+  myFirstParameter = theFirst;
-  DomainSup=DSup;
+  myLastParameter = theLast;
 }

src/IntAna/IntAna_Curve.hxx

@@ -21,16 +21,9 @@
 #include <Standard_DefineAlloc.hxx>
 #include <Standard_Handle.hxx>
 
-#include <Standard_Real.hxx>
-#include <Standard_Boolean.hxx>
 #include <GeomAbs_SurfaceType.hxx>
 #include <gp_Ax3.hxx>
-class Standard_DomainError;
+#include <TColStd_ListOfReal.hxx>
-class gp_Cylinder;
-class gp_Cone;
-class gp_Pnt;
-class gp_Vec;
-
 
 //! Definition of a parametric Curve which is the result
 //! of the intersection between two quadrics.
@@ -57,7 +50,7 @@ public:
   Standard_EXPORT Standard_Boolean IsOpen() const;
   
   //! Returns the paramatric domain of the curve.
-  Standard_EXPORT void Domain (Standard_Real& Theta1, Standard_Real& Theta2) const;
+  Standard_EXPORT void Domain(Standard_Real& theFirst, Standard_Real& theLast) const;
   
   //! Returns TRUE if the function is constant.
   Standard_EXPORT Standard_Boolean IsConstant() const;
@@ -77,11 +70,11 @@ public:
   
   //! Tries to find the parameter of the point P on the curve.
   //! If the method returns False, the "projection" is
-  //! impossible, and the value of Para is not significant.
+  //! impossible.
-  //! If the method returns True, Para is the parameter of the
+  //! If the method returns True at least one parameter has been found.
-  //! nearest intersection between the curve and the iso-theta
+  //! theParams is always sorted in ascending order.
-  //! containing P.
+  Standard_EXPORT void FindParameter(const gp_Pnt& P,
-  Standard_EXPORT Standard_Boolean FindParameter (const gp_Pnt& P, Standard_Real& Para) const;
+                                     TColStd_ListOfReal& theParams) const;
   
   //! If flag is True, the Curve is not defined at the
   //! first parameter of its domain.
@@ -91,10 +84,8 @@ public:
   //! first parameter of its domain.
   Standard_EXPORT void SetIsLastOpen (const Standard_Boolean Flag);
   
-  //! Protected function.
+  //! Trims this curve
-  Standard_EXPORT void InternalUVValue (const Standard_Real Param, Standard_Real& U, Standard_Real& V, Standard_Real& A, Standard_Real& B, Standard_Real& C, Standard_Real& Co, Standard_Real& Si, Standard_Real& Di) const;
+  Standard_EXPORT void SetDomain(const Standard_Real theFirst, const Standard_Real theLast);
-  
-  Standard_EXPORT void SetDomain (const Standard_Real Theta1, const Standard_Real Theta2);
 
 
 
@@ -105,6 +96,8 @@ protected:
   //! Protected function.
   Standard_EXPORT gp_Pnt InternalValue (const Standard_Real Theta1, const Standard_Real Theta2) const;
 
+  //! Protected function.
+  Standard_EXPORT void InternalUVValue (const Standard_Real Param, Standard_Real& U, Standard_Real& V, Standard_Real& A, Standard_Real& B, Standard_Real& C, Standard_Real& Co, Standard_Real& Si, Standard_Real& Di) const;
 
 
 
@@ -133,8 +126,9 @@ private:
   Standard_Boolean TwoCurves;
   Standard_Boolean TakeZPositive;
   Standard_Real Tolerance;
-  Standard_Real DomainInf;
+
-  Standard_Real DomainSup;
+  //! Internal fields defining the default domain
+  Standard_Real DomainInf, DomainSup;
   Standard_Boolean RestrictedInf;
   Standard_Boolean RestrictedSup;
   Standard_Boolean firstbounded;
@@ -144,6 +138,9 @@ private:
   Standard_Real Angle;
   gp_Ax3 Ax3;
 
+  //! Trim boundaries
+  Standard_Real myFirstParameter, myLastParameter;
+
 
 };
 

src/IntAna/IntAna_IntQuadQuad.cxx

@@ -28,6 +28,7 @@
 //==                           C Y L I N D R E   Q U A D R I Q U E
 //======================================================================
 
+#include <ElSLib.hxx>
 #include <gp_Ax2.hxx>
 #include <gp_Ax3.hxx>
 #include <gp_Cone.hxx>
@@ -41,6 +42,75 @@
 #include <Standard_OutOfRange.hxx>
 #include <StdFail_NotDone.hxx>
 
+//=======================================================================
+//function : AddSpecialPoints
+//purpose  : Sometimes the boundaries theTheta1 and theTheta2 are
+//            computed with some inaccuracy. At that, some special points
+//            (cone apex or sphere pole(s)), which are true intersection
+//            points lie out of the domain [theTheta1, theTheta2] of the ALine.
+//           This function corrects these boundaries to make them be included 
+//            in the domain of the ALine.
+//           Parameters Theta1 and Theta2 must be initialized
+//            before calling this function.
+//=======================================================================
+template <class gpSmth>
+static void AddSpecialPoints(const IntAna_Quadric& theQuad,
+                             const gpSmth& theGpObj,
+                             Standard_Real& theTheta1,
+                             Standard_Real& theTheta2)
+{
+  const Standard_Real aPeriod = M_PI + M_PI;
+  const NCollection_List<gp_Pnt> &aLSP = theQuad.SpecialPoints();
+
+  if (aLSP.IsEmpty())
+    return;
+
+  Standard_Real aU = 0.0, aV = 0.0;
+  Standard_Real aMaxDelta = 0.0;
+  for (NCollection_List<gp_Pnt>::Iterator anItr(aLSP); anItr.More(); anItr.Next())
+  {
+    const gp_Pnt &aPt = anItr.Value();
+    ElSLib::Parameters(theGpObj, aPt, aU, aV);
+    const gp_Pnt aPProj(ElSLib::Value(aU, aV, theGpObj));
+
+    if (aPt.SquareDistance(aPProj) > Precision::SquareConfusion())
+    {
+      // aPt is not an intersection point
+      continue;
+    }
+
+    Standard_Real aDelta1 = Min(aU - theTheta1, 0.0),
+                  aDelta2 = Max(aU - theTheta2, 0.0);
+
+    if (aDelta1 < -M_PI)
+    {
+      // Must be aDelta1 = Min(aU - theTheta1 + aPeriod, 0.0).
+      // But aU - theTheta1 + aPeriod >= 0 always.
+      aDelta1 = 0.0;
+    }
+
+    if (aDelta2 > M_PI)
+    {
+      // Must be aDelta2 = Max(aU - theTheta2 - aPeriod, 0.0).
+      // But aU - theTheta2 - aPeriod <= 0 always.
+      aDelta2 = 0.0;
+    }
+
+    const Standard_Real aDelta = Max(-aDelta1, aDelta2);
+    aMaxDelta = Max(aMaxDelta, aDelta);
+  }
+
+  if(aMaxDelta != 0.0)
+  {
+    theTheta1 -= aMaxDelta;
+    theTheta2 += aMaxDelta;
+    if ((theTheta2 - theTheta1) > aPeriod)
+    {
+      theTheta2 = theTheta1 + aPeriod;
+    }
+  }
+}
+
 //=======================================================================
 //class : TrigonometricRoots
 //purpose: Classe Interne (Donne des racines classees d un polynome trigo)
@@ -486,13 +556,15 @@ void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
 		//
 		qwet=MTF.Value(autrepar);
 		if(qwet>=0.) { 
+                  Standard_Real aParam = Theta1 + PIpPI;
+                  AddSpecialPoints(Quad, Cyl, Theta1, aParam);
 		  TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
-							   myEpsilon,Theta1,Theta1+PIpPI,
+							   myEpsilon,Theta1,aParam,
 							   UN_SEUL_Z_PAR_THETA,
 							   Z_POSITIF);
 		  NbCurves++;
 		  TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
-							   myEpsilon,Theta1,Theta1+PIpPI,
+							   myEpsilon,Theta1,aParam,
 							   UN_SEUL_Z_PAR_THETA,
 							   Z_NEGATIF);
 		  NbCurves++;
@@ -534,6 +606,7 @@ void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
 		//ft
 		if((Theta3-Theta2)<5.e-8) { 
 		//
+                  AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
 		  TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 							   myEpsilon,Theta1,Theta2,
 							   UN_SEUL_Z_PAR_THETA,
@@ -546,6 +619,7 @@ void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
 		  NbCurves++;
 		}
 		else { 
+                  AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
 		  TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 							   myEpsilon,Theta1,Theta2,
 							   DEUX_Z_PAR_THETA,

src/IntAna/IntAna_Quadric.cxx

@@ -28,6 +28,7 @@
 #include <gp_Trsf.hxx>
 #include <gp_Vec.hxx>
 #include <IntAna_Quadric.hxx>
+#include <ElSLib.hxx>
 
 //----------------------------------------------------------------------
 //-- 
@@ -85,22 +86,27 @@ IntAna_Quadric::IntAna_Quadric(const gp_Cylinder& Cyl) {
 //----------------------------------------------------------------------
 //-- Cone -----> Quadric
 //----------------------------------------------------------------------
-IntAna_Quadric::IntAna_Quadric(const gp_Cone& Cone) {
+IntAna_Quadric::IntAna_Quadric(const gp_Cone& Cone)
-  Cone.Coefficients(CXX,CYY,CZZ,CXY,CXZ,CYZ,CX,CY,CZ,CCte);  
+{
+  SetQuadric(Cone);
 }
 
 void IntAna_Quadric::SetQuadric(const gp_Cone& Cone) {
   Cone.Coefficients(CXX,CYY,CZZ,CXY,CXZ,CYZ,CX,CY,CZ,CCte);  
+  const Standard_Real aVParam = -Cone.RefRadius() / Sin(Cone.SemiAngle());
+  mySpecialPoints.Append(ElSLib::Value(0.0, aVParam, Cone));
 }
 //----------------------------------------------------------------------
 //-- Sphere -----> Quadric
 //----------------------------------------------------------------------
 void IntAna_Quadric::SetQuadric(const gp_Sphere& Sph) {
   Sph.Coefficients(CXX,CYY,CZZ,CXY,CXZ,CYZ,CX,CY,CZ,CCte);
+  mySpecialPoints.Append(ElSLib::Value(0.0, -M_PI_2, Sph));
+  mySpecialPoints.Append(ElSLib::Value(0.0, M_PI_2, Sph));
 }
 
 IntAna_Quadric::IntAna_Quadric(const gp_Sphere& Sph) {
-  Sph.Coefficients(CXX,CYY,CZZ,CXY,CXZ,CYZ,CX,CY,CZ,CCte);
+  SetQuadric(Sph);
 }
 //----------------------------------------------------------------------
 //-- Returns the Coefficients of the Quadric

src/IntAna/IntAna_Quadric.hxx

@@ -17,17 +17,8 @@
 #ifndef _IntAna_Quadric_HeaderFile
 #define _IntAna_Quadric_HeaderFile
 
-#include <Standard.hxx>
 #include <Standard_DefineAlloc.hxx>
-#include <Standard_Handle.hxx>
+#include <NCollection_List.hxx>
-
-#include <Standard_Real.hxx>
-class gp_Pln;
-class gp_Sphere;
-class gp_Cylinder;
-class gp_Cone;
-class gp_Ax3;
-
 
 //! This class provides a description of Quadrics by their
 //! Coefficients in natural coordinate system.
@@ -78,7 +69,11 @@ public:
   //! in the local coordinates system defined by Axis
   Standard_EXPORT void NewCoefficients (Standard_Real& xCXX, Standard_Real& xCYY, Standard_Real& xCZZ, Standard_Real& xCXY, Standard_Real& xCXZ, Standard_Real& xCYZ, Standard_Real& xCX, Standard_Real& xCY, Standard_Real& xCZ, Standard_Real& xCCte, const gp_Ax3& Axis) const;
 
-
+  //! Returns the list of special points (with singularities)
+  const NCollection_List<gp_Pnt>& SpecialPoints() const
+  {
+    return mySpecialPoints;
+  }
 
 
 protected:
@@ -101,7 +96,7 @@ private:
   Standard_Real CY;
   Standard_Real CZ;
   Standard_Real CCte;
-
+  NCollection_List<gp_Pnt> mySpecialPoints;
 
 };
 

src/IntPatch/IntPatch_ALine.hxx

@@ -17,31 +17,21 @@
 #ifndef _IntPatch_ALine_HeaderFile
 #define _IntPatch_ALine_HeaderFile
 
-#include <Standard.hxx>
+#include <Standard_Handle.hxx>
-#include <Standard_Type.hxx>
 
 #include <IntAna_Curve.hxx>
-#include <Standard_Boolean.hxx>
-#include <Standard_Integer.hxx>
-#include <IntPatch_SequenceOfPoint.hxx>
 #include <IntPatch_Line.hxx>
-#include <IntSurf_TypeTrans.hxx>
+#include <IntPatch_SequenceOfPoint.hxx>
-#include <IntSurf_Situation.hxx>
+#include <TColStd_ListOfReal.hxx>
-#include <Standard_Real.hxx>
+
-#include <gp_Pnt.hxx>
-class Standard_DomainError;
-class Standard_OutOfRange;
 class IntAna_Curve;
 class IntPatch_Point;
-class gp_Pnt;
-class gp_Vec;
-
-
 class IntPatch_ALine;
+
 DEFINE_STANDARD_HANDLE(IntPatch_ALine, IntPatch_Line)
 
 //! Implementation of an intersection line described by a
-//! parametrised curve.
+//! parametrized curve.
 class IntPatch_ALine : public IntPatch_Line
 {
 
@@ -97,13 +87,13 @@ public:
   //! intersection.
     Standard_Boolean D1 (const Standard_Real U, gp_Pnt& P, gp_Vec& Du);
   
-  //! Tries to find the parameter of the point P on the curve.
+  //! Tries to find the parameters of the point P on the curve.
   //! If the method returns False, the "projection" is
-  //! impossible, and the value of Para is not significant.
+  //! impossible.
-  //! If the method returns True, Para is the parameter of the
+  //! If the method returns True at least one parameter has been found.
-  //! nearest intersection between the curve and the iso-theta
+  //! theParams is always sorted in ascending order.
-  //! containing P.
+  void FindParameter(const gp_Pnt& P,
-    Standard_Boolean FindParameter (const gp_Pnt& P, Standard_Real& Para) const;
+                     TColStd_ListOfReal& theParams) const;
   
   //! Returns True if the line has a known First point.
   //! This point is given by the method FirstPoint().
@@ -126,6 +116,12 @@ public:
   //! Returns the vertex of range Index on the line.
     const IntPatch_Point& Vertex (const Standard_Integer Index) const;
   
+  //! Allows modifying the vertex with index theIndex on the line.
+  IntPatch_Point& ChangeVertex(const Standard_Integer theIndex)
+  {
+    return svtx.ChangeValue(theIndex);
+  }
+
   //! Set the parameters of all the vertex on the line.
   //! if a vertex is already in the line,
   //! its parameter is modified

src/IntPatch/IntPatch_ALine.lxx

@@ -62,9 +62,10 @@ inline Standard_Boolean IntPatch_ALine::D1(const Standard_Real U, gp_Pnt& P, gp_
   return curv.D1u(U,P,Du); // D1u leve l exception DomainError
 }
 
-inline Standard_Boolean IntPatch_ALine::FindParameter(const gp_Pnt& P, Standard_Real& Para) const
+inline void IntPatch_ALine::FindParameter(const gp_Pnt& theP,
+                                          TColStd_ListOfReal& theParams) const
 {
-  return curv.FindParameter(P,Para);
+  curv.FindParameter(theP, theParams);
 }
 
 inline Standard_Boolean IntPatch_ALine::HasFirstPoint () const

src/IntPatch/IntPatch_ALineToWLine.cxx

@@ -62,21 +62,20 @@ static void AddVertexPoint(Handle(IntSurf_LineOn2S)& theLine,
 //function : IsPoleOrSeam
 //purpose  : Processes theVertex depending on its type
 //            (pole/apex/point on boundary etc.) and adds it in theLine.
+//           thePIsoRef is the reference point using in case when the
+//            value of correspond parameter cannot be precise.
 //           theSingularSurfaceID contains the ID of surface with
 //            special point (0 - none, 1 - theS1, 2 - theS2)
 //=======================================================================
 static IntPatch_SpecPntType IsPoleOrSeam(const Handle(Adaptor3d_HSurface)& theS1,
                                          const Handle(Adaptor3d_HSurface)& theS2,
+                                         const IntSurf_PntOn2S& thePIsoRef,
                                          Handle(IntSurf_LineOn2S)& theLine,
                                          IntPatch_Point &theVertex,
-                                         const Standard_Real* const theArrPeriods,
+                                         const Standard_Real theArrPeriods[4],
                                          const Standard_Real theTol3d,
                                          Standard_Integer& theSingularSurfaceID)
 {
-  const Standard_Integer aNbPnts = theLine->NbPoints();
-  if(aNbPnts == 0)
-    return IntPatch_SPntNone;
-
   theSingularSurfaceID = 0;
 
   for(Standard_Integer i = 0; i < 2; i++)
@@ -94,8 +93,8 @@ static IntPatch_SpecPntType IsPoleOrSeam(const Handle(Adaptor3d_HSurface)& theS1
       {
         if(IntPatch_SpecialPoints::
               AddSingularPole((isReversed? theS2 : theS1), (isReversed? theS1 : theS2),
-                               theLine->Value(aNbPnts), theTol3d, theVertex,
+                              thePIsoRef, theVertex, anApexPoint,
-                               anApexPoint, isReversed, Standard_True))
+                              isReversed, Standard_True))
         {
           anAddedPType = IntPatch_SPntPole;
           break;
@@ -107,8 +106,8 @@ static IntPatch_SpecPntType IsPoleOrSeam(const Handle(Adaptor3d_HSurface)& theS1
       {
         if(IntPatch_SpecialPoints::
             AddCrossUVIsoPoint((isReversed? theS2 : theS1), (isReversed? theS1 : theS2),
-                                          theLine->Value(aNbPnts), theTol3d,
+                                  thePIsoRef, theTol3d,
-                                          anApexPoint, isReversed))
+                                  anApexPoint, isReversed))
         {
           anAddedPType = IntPatch_SPntSeamUV;
           break;
@@ -208,7 +207,7 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 //function : SetTol3D
 //purpose  : 
 //=======================================================================
-  void IntPatch_ALineToWLine::SetTol3D(const Standard_Real aTol)
+void IntPatch_ALineToWLine::SetTol3D(const Standard_Real aTol)
 {
   myTol3D = aTol;
 }
@@ -216,7 +215,7 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 //function : Tol3D
 //purpose  : 
 //=======================================================================
-  Standard_Real IntPatch_ALineToWLine::Tol3D()const
+Standard_Real IntPatch_ALineToWLine::Tol3D()const
 {
   return myTol3D;
 }
@@ -224,7 +223,7 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 //function : SetTolTransition
 //purpose  : 
 //=======================================================================
-  void IntPatch_ALineToWLine::SetTolTransition(const Standard_Real aTol)
+void IntPatch_ALineToWLine::SetTolTransition(const Standard_Real aTol)
 {
   myTolTransition = aTol;
 }
@@ -232,7 +231,7 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 //function : TolTransition
 //purpose  : 
 //=======================================================================
-  Standard_Real IntPatch_ALineToWLine::TolTransition()const
+Standard_Real IntPatch_ALineToWLine::TolTransition()const
 {
   return myTolTransition;
 }
@@ -240,7 +239,7 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 //function : SetTolOpenDomain
 //purpose  : 
 //=======================================================================
-  void IntPatch_ALineToWLine::SetTolOpenDomain(const Standard_Real aTol)
+void IntPatch_ALineToWLine::SetTolOpenDomain(const Standard_Real aTol)
 {
   myTolOpenDomain = aTol;
 }
@@ -252,6 +251,48 @@ IntPatch_ALineToWLine::IntPatch_ALineToWLine(const Handle(Adaptor3d_HSurface)& t
 {
   return myTolOpenDomain;
 }
+
+//=======================================================================
+//function : GetSectionRadius
+//purpose  : 
+//=======================================================================
+Standard_Real IntPatch_ALineToWLine::GetSectionRadius(const gp_Pnt& thePnt3d) const
+{
+  Standard_Real aRetVal = RealLast();
+  for (Standard_Integer i = 0; i < 2; i++)
+  {
+    const IntSurf_Quadric& aQuad = i ? myQuad2 : myQuad1;
+    if (aQuad.TypeQuadric() == GeomAbs_Cone)
+    {
+      const gp_Cone aCone = aQuad.Cone();
+      const gp_XYZ aRVec = thePnt3d.XYZ() - aCone.Apex().XYZ();
+      const gp_XYZ &aDir = aCone.Axis().Direction().XYZ();
+
+      aRetVal = Min(aRetVal, Abs(aRVec.Dot(aDir)*Tan(aCone.SemiAngle())));
+    }
+    else if (aQuad.TypeQuadric() == GeomAbs_Sphere)
+    {
+      const gp_Sphere aSphere = aQuad.Sphere();
+      const gp_XYZ aRVec = thePnt3d.XYZ() - aSphere.Location().XYZ();
+      const gp_XYZ &aDir = aSphere.Position().Direction().XYZ();
+      const Standard_Real aR = aSphere.Radius();
+      const Standard_Real aD = aRVec.Dot(aDir);
+      const Standard_Real aDelta = aR*aR - aD*aD;
+      if (aDelta <= 0.0)
+      {
+        aRetVal = 0.0;
+        break;
+      }
+      else
+      {
+        aRetVal = Min(aRetVal, Sqrt(aDelta));
+      }
+    }
+  }
+
+  return aRetVal;
+}
+
 //=======================================================================
 //function : MakeWLine
 //purpose  : 
@@ -282,10 +323,81 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
                                       IntPatch_SequenceOfLine& theLines) const 
 {
   const Standard_Integer aNbVert = theALine->NbVertex();
-  if (!aNbVert) {
+  if (aNbVert == 0)
+  {
     return;
   }
+
+#if 0
+  //To draw ALine as a wire DRAW-object use the following code.
+  {
+    static int zzz = 0;
+    zzz++;
+
+    bool flShow = /*(zzz == 1)*/false;
+
+    if (flShow)
+    {
+      std::cout << " +++ DUMP ALine (begin) +++++" << std::endl;
+      Standard_Integer aI = 0;
+      const Standard_Real aStep = (theLPar - theFPar) / 9999.0;
+      for (Standard_Real aPrm = theFPar; aPrm < theLPar; aPrm += aStep)
+      {
+        const gp_Pnt aPP(theALine->Value(aPrm));
+        std::cout << "vertex v" << ++aI << " " << aPP.X() << " " << aPP.Y() << " " << aPP.Z() << std::endl;
+      }
+
+      gp_Pnt aPP(theALine->Value(theLPar));
+      std::cout << "vertex v" << ++aI << " " << aPP.X() << " " << aPP.Y() << " " << aPP.Z() << std::endl;
+      std::cout << " --- DUMP ALine (end) -----" << std::endl;
+    }
+  }
+
+  //Copy all output information and apply it as a TCL-code in DRAW.
+
+  //After that, use TCL-script below:
+
+  /* ********************************* Script (begin)
+  shape ww w
+  copy v1 vprev
+  for {set i 2} {$i <= 10000} {incr i} {
+    distmini dd vprev v$i;
+
+    if { [dval dd_val] > 1.0e-7} {
+      edge ee vprev v$i;
+      add ee ww;
+      copy v$i vprev;
+    }
+  }
+  ********************************** Script (end) */
+#endif
+
+  //The same points can be marked by different vertices.
+  //The code below unifies tolerances of all vertices
+  //marking the same points.
+  for (Standard_Integer i = 1; i < aNbVert; i++)
+  {
+    IntPatch_Point &aCurVert = theALine->ChangeVertex(i);
+    const IntSurf_PntOn2S &aCurrPt = aCurVert.PntOn2S();
+    const Standard_Real aCurToler = aCurVert.Tolerance();
+    for (Standard_Integer j = i + 1; j <= aNbVert; j++)
+    {
+      IntPatch_Point &aVert = theALine->ChangeVertex(j);
+      const IntSurf_PntOn2S &aNewPt = aVert.PntOn2S();
+      const Standard_Real aToler = aVert.Tolerance();
+
+      const Standard_Real aSumTol = aCurToler + aToler;
+      if (aCurrPt.IsSame(aNewPt, aSumTol))
+      {
+        aCurVert.SetTolerance(aSumTol);
+        aVert.SetTolerance(aSumTol);
+      }
+    }
+  }
+
   const Standard_Real aTol = 2.0*myTol3D+Precision::Confusion();
+  const Standard_Real aPrmTol = Max(1.0e-4*(theLPar - theFPar), Precision::PConfusion());
+
   IntPatch_SpecPntType aPrePointExist = IntPatch_SPntNone;
   
   NCollection_Array1<Standard_Real> aVertexParams(1, aNbVert);
@@ -299,7 +411,8 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
   
   for(Standard_Integer i = aVertexParams.Lower(); i <= aVertexParams.Upper(); i++)
   {
-    const Standard_Real aPar = theALine->Vertex(i).ParameterOnLine();
+    const IntPatch_Point& aVert = theALine->Vertex(i);
+    const Standard_Real aPar = aVert.ParameterOnLine();
     aVertexParams(i) = aPar;
     hasVertexBeenChecked(i) = Standard_False;
   }
@@ -344,11 +457,24 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
         }
       }
 
+      Standard_Boolean isPointValid = Standard_False;
       Standard_Real aTgMagn = 0.0;
       {
         gp_Pnt aPnt3d;
         gp_Vec aTg;
         theALine->D1(aParameter, aPnt3d, aTg);
+        if (GetSectionRadius(aPnt3d) < 5.0e-6)
+        {
+          // We cannot compute 2D-parameters of
+          // aPOn2S correctly.
+
+          isPointValid = Standard_False;
+        }
+        else
+        {
+          isPointValid = Standard_True;
+        }
+        
         aTgMagn = aTg.Magnitude();
         Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
         myQuad1.Parameters(aPnt3d, u1, v1);
@@ -372,13 +498,25 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
         if (aPrePointExist == IntPatch_SPntPole)
         {
           Standard_Real aPrt = 0.5*(aPrevParam + theLPar);
-          for (Standard_Integer i = aVertexParams.Lower(); i <= aVertexParams.Upper(); i++)
+          for (Standard_Integer i = aVertexParams.Lower();
+               i <= aVertexParams.Upper(); i++)
           {
             const Standard_Real aParam = aVertexParams(i);
 
             if (aParam <= aPrevParam)
               continue;
 
+            if ((aParam - aPrevParam) < aPrmTol)
+            {
+              const gp_Pnt aPnt3d(theALine->Value(aParam));
+              if (aPOn2S.Value().SquareDistance(aPnt3d) < Precision::SquareConfusion())
+              {
+                // i-th vertex is the same as a Pole/Apex.
+                // So, it should be ignored.
+                continue;
+              }
+            }
+
             aPrt = 0.5*(aParam + aPrevParam);
             break;
           }
@@ -426,10 +564,17 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
         if(hasVertexBeenChecked(i))
           continue;
 
+        const IntPatch_Point &aVP = theALine->Vertex(i);
         const Standard_Real aParam = aVertexParams(i);
         if( ((aPrevParam < aParam) && (aParam <= aParameter)) ||
-            ((aPrevParam == aParameter) && (aParam == aParameter)))
+            ((aPrevParam == aParameter) && (aParam == aParameter))||
+            (aPOn2S.IsSame(aVP.PntOn2S(), aVP.Tolerance()) && 
+                    (Abs(aVP.ParameterOnLine() - aParameter) < aPrmTol)))
         {
+          //We have either jumped over the vertex or "fell" on the vertex.
+          //However, ALine can be self-interfered. Therefore, we need to check
+          //vertex parameter and 3D-distance together.
+
           aVertexNumber = i;
           break;
         }
@@ -439,10 +584,14 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
       
       if(aVertexNumber < 0)
       {
-        StepComputing(theALine, aPOn2S, theLPar, aParameter, aTgMagn,
+        if (isPointValid)
-                              aStepMin, aStepMax, myTol3D, aStep);
+        {
-        AddPointIntoLine(aLinOn2S, anArrPeriods, aPOn2S);
+          StepComputing(theALine, aPOn2S, theLPar, aParameter, aTgMagn,
-        aPrevLPoint = aPOn2S;
+                        aStepMin, aStepMax, myTol3D, aStep);
+          AddPointIntoLine(aLinOn2S, anArrPeriods, aPOn2S);
+          aPrevLPoint = aPOn2S;
+        }
+
         continue;
       }
 
@@ -453,7 +602,33 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
       // IsPoleOrSeam inserts new point in aLinOn2S if aVtx respects
       //to some special point. Otherwise, aLinOn2S is not changed.
 
-      aPrePointExist = IsPoleOrSeam(myS1, myS2, aLinOn2S, aVtx,
+      // Find a point for reference parameter. It will be used
+      // if real parameter value cannot be precise (see comment to 
+      // IsPoleOrSeam(...) function). 
+      IntSurf_PntOn2S aPrefIso = aVtx.PntOn2S();
+      if (aLinOn2S->NbPoints() < 1)
+      {
+        for (Standard_Integer i = aVertexNumber + 1; i <= aVertexParams.Upper(); i++)
+        {
+          const Standard_Real aParam = aVertexParams(i);
+          if ((aParam - aVertexParams(aVertexNumber)) > Precision::PConfusion())
+          {
+            const Standard_Real aPrm = 0.5*(aParam + aVertexParams(aVertexNumber));
+            const gp_Pnt aPnt3d(theALine->Value(aPrm));
+            Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
+            myQuad1.Parameters(aPnt3d, u1, v1);
+            myQuad2.Parameters(aPnt3d, u2, v2);
+            aPrefIso.SetValue(aPnt3d, u1, v1, u2, v2);
+            break;
+          }
+        }
+      }
+      else
+      {
+        aPrefIso = aLinOn2S->Value(aLinOn2S->NbPoints());
+      }
+
+      aPrePointExist = IsPoleOrSeam(myS1, myS2, aPrefIso, aLinOn2S, aVtx,
                                 anArrPeriods, aTol, aSingularSurfaceID);
 
       const Standard_Real aCurVertParam = aVtx.ParameterOnLine();
@@ -463,6 +638,12 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
       }
       else
       {
+        if (!isPointValid)
+        {
+          //Take a farther point of ALine (with greater parameter)
+          continue;
+        }
+
         if(aVtx.Tolerance() > aTol)
         {
           aVtx.SetValue(aPOn2S);
@@ -507,7 +688,7 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
         }
       }
 
-      if(aPrePointExist != IntPatch_SPntNone)
+      if ((aPrePointExist != IntPatch_SPntNone) && (aLinOn2S->NbPoints() > 1))
         break;
     }//for(; !isLast; aParameter += aStep)
 
@@ -586,6 +767,9 @@ void IntPatch_ALineToWLine::MakeWLine(const Handle(IntPatch_ALine)& theALine,
     if (aWLine->NbPnts() > 1)
     {
       aWLine->EnablePurging(Standard_False);
+#ifdef INTPATCH_ALINETOWLINE_DEBUG
+      aWLine->Dump(0);
+#endif
       theLines.Append(aWLine);
     }
   }//while(aParameter < theLPar)
@@ -653,6 +837,19 @@ Standard_Boolean IntPatch_ALineToWLine::
   const Standard_Real aR = IntPatch_PointLine::
                             CurvatureRadiusOfIntersLine(myS1, myS2, thePOn2S);
 
+#if 0
+  {
+    static int zzz = 0;
+    zzz++;
+    std::cout << "*** R" << zzz << " (begin)" << std::endl;
+    Standard_Real aU1, aV1, aU2, aV2;
+    thePOn2S.Parameters(aU1, aV1, aU2, aV2);
+    std::cout << "Prms: " << aU1 << ", " << aV1 << ", " << aU2 << ", " << aV2 << std::endl;
+    std::cout << "Radius = " << aR << std::endl;
+    std::cout << "*** R" << zzz << " (end)" << std::endl;
+  }
+#endif
+
   if(aR < 0.0)
   {
     return Standard_False;

src/IntPatch/IntPatch_ALineToWLine.hxx

@@ -80,11 +80,18 @@ protected:
                                                  const Standard_Real theMaxDeflection,
                                                  Standard_Real& theStep) const;
 
-
+  //! Compares distances from theMidPt to every quadrics with theMaxDeflection
+  //! (maximal distance of two ones is taken into account).
+  //! Returns the result of this comparison: -1 - small distance, +1 - big distance,
+  //! 0 - Dist == theMaxDeflection. Comparisons are done with internal tolerances.
   Standard_EXPORT Standard_Integer CheckDeflection(const gp_XYZ& theMidPt,
                                                    const Standard_Real theMaxDeflection) const;
 
-
+  //! Returns radius of a circle obtained by intersection the quadric with a plane
+  //! goes through thePnt3d perpendicular to the quadric axis. This radius is computed
+  //! for both quadrics and minimal value is returned.
+  //! This check is made for cone and sphere only.
+  Standard_EXPORT Standard_Real GetSectionRadius(const gp_Pnt& thePnt3d) const;
 
 private:
 

src/IntPatch/IntPatch_ImpImpIntersection_0.gxx

@@ -42,24 +42,26 @@ static void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 			    const Standard_Boolean,
 			    const Standard_Real);
 
-static Standard_Boolean MultiplePoint (const IntPatch_SequenceOfPathPointOfTheSOnBounds&,
+static Standard_Boolean MultiplePoint(const IntPatch_SequenceOfPathPointOfTheSOnBounds& listpnt,
-				       const Handle(Adaptor3d_TopolTool)&,
+                                      const Handle(Adaptor3d_TopolTool)& Domain,
-				       const IntSurf_Quadric&,
+                                      const IntSurf_Quadric& QuadSurf,
-		                       const gp_Vec&,
+                                      const gp_Vec&    Normale,
-		                       const IntPatch_SequenceOfLine&,
+                                      const IntPatch_SequenceOfLine& slin,
-                                       TColStd_Array1OfInteger&,
+                                      TColStd_Array1OfInteger& Done,
-                                       TColStd_Array1OfInteger&,
+                                      TColStd_Array1OfInteger& UsedLine,
-                                       const Standard_Integer,
+                                      const Standard_Integer Index,
-		                       const Standard_Boolean);
+                                      const Standard_Boolean OnFirst,
+                                      const Standard_Real theToler);
 
-static Standard_Boolean PointOnSecondDom (const IntPatch_SequenceOfPathPointOfTheSOnBounds&,
+static Standard_Boolean PointOnSecondDom(const IntPatch_SequenceOfPathPointOfTheSOnBounds& listpnt,
-					  const Handle(Adaptor3d_TopolTool)&,
+                                         const Handle(Adaptor3d_TopolTool)& Domain,
-					  const IntSurf_Quadric&,
+                                         const IntSurf_Quadric& QuadSurf,
-		                          const gp_Vec&,
+                                         const gp_Vec&    Normale,
-		                          const gp_Vec&,
+                                         const gp_Vec&    Vtgint,
-		                          const Handle(IntPatch_Line)&,
+                                         const Handle(IntPatch_Line)& lin,
-                                          TColStd_Array1OfInteger&,
+                                         TColStd_Array1OfInteger& Done,
-                                          const Standard_Integer);
+                                         const Standard_Integer Index,
+                                         const Standard_Real theToler);
 
 static Standard_Boolean SingleLine (const gp_Pnt&,
 				    const Handle(IntPatch_Line)&,
@@ -68,17 +70,19 @@ static Standard_Boolean SingleLine (const gp_Pnt&,
 				    gp_Vec&);
 
 
-static Standard_Boolean FindLine (gp_Pnt&,
+static Standard_Boolean FindLine(gp_Pnt& Psurf,
-				  const IntPatch_SequenceOfLine&,
+                                 const IntPatch_SequenceOfLine& slin,
-				  const Standard_Real,
+                                 const Standard_Real Tol,
-				  Standard_Real&,
+                                 TColStd_ListOfReal& theLParams,
-				  gp_Vec&,
+                                 gp_Vec& Vtgtint,
-				  Standard_Integer&,
+                                 Standard_Integer& theLineIdx,
-				  Standard_Integer,
+                                 Standard_Integer OnlyThisLine,
-				  const Handle(Adaptor2d_HCurve2d)&,
+                                 const Handle(Adaptor2d_HCurve2d)& thearc,
-				  Standard_Real&,
+                                 Standard_Real& theparameteronarc,
-				  gp_Pnt& pointonarc,
+                                 gp_Pnt& thepointonarc,
-				  const IntSurf_Quadric&);
+                                 const IntSurf_Quadric& QuadSurf1,
+                                 const IntSurf_Quadric& QuadSurf2,
+                                 Standard_Real& theOutputToler);
 
 static void ProcessSegments (const IntPatch_SequenceOfSegmentOfTheSOnBounds&,
 			     IntPatch_SequenceOfLine&,
@@ -102,8 +106,7 @@ Standard_Boolean IntersectionWithAnArc(gp_Pnt& PSurf,
 				       gp_Pnt& thepointonarc,
 				       const IntSurf_Quadric& QuadSurf,
 				       const Standard_Real u0alin,
-				       const Standard_Real u1alin,
+				       const Standard_Real u1alin) { 
-				       Standard_Real& actualdist) { 
   Standard_Real dtheta,theta;
 #ifdef OCCT_DEBUG
   //Standard_Real u,v,A,B,C,cost,sint,sign;
@@ -221,7 +224,6 @@ Standard_Boolean IntersectionWithAnArc(gp_Pnt& PSurf,
       thepointonarc = alin->Value(para);
       cpasok=Standard_False;
 //--      printf("\nt:%d",nbiter);
-      actualdist = bestdist;
       return(Standard_True);
     }
     else { 
@@ -249,7 +251,6 @@ Standard_Boolean IntersectionWithAnArc(gp_Pnt& PSurf,
     _theparameteronarc=bestpara;
     thepointonarc = alin->Value(para);
 //--     printf("\nT:%d",nbiter);
-    actualdist=bestdist;
     return(Standard_True);
   }
 //--   printf("\nF:%d",nbiter);
@@ -383,7 +384,7 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
   //
   Standard_Integer i,k;
   Standard_Integer linenumber;
-  Standard_Real paraint = 0.,currentparameter,tolerance;
+  Standard_Real currentparameter,tolerance;
   Standard_Real U1,V1,U2,V2;
   Standard_Boolean goon;
   
@@ -429,10 +430,10 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 	Vtgrst.SetLinearForm(d2d.X(),d1u,d2d.Y(),d1v);
 #endif      
         goon = MultiplePoint(listpnt,Domain,QuadSurf,Normale,slin,Done, UsedLine,
-			     i,OnFirst);
+                             i, OnFirst, Tolarc);
       }
       if (goon) {
-	Standard_Boolean linefound;
+        Standard_Boolean linefound = Standard_False;
 	
 	for(Standard_Integer indiceline = 1; indiceline <=slin.Length(); indiceline++) { 
 	  if( UsedLine(indiceline) != 0 )
@@ -470,8 +471,11 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 	  //  Modified by skv - Thu Jan 15 15:57:15 2004 OCC4455 End
 	  gp_Pnt pointonarc;
 	  Vtgint.SetCoord(0,0,0);
-	  linefound = FindLine(Psurf,slin,tolerance,paraint,Vtgint,linenumber,indiceline,
+          Standard_Real aVertTol = Tolarc;
-			       currentarc,currentparameter,pointonarc,QuadSurf);  
+          TColStd_ListOfReal aLParams;
+          linefound = FindLine(Psurf, slin, tolerance, aLParams, Vtgint, linenumber,
+                               indiceline, currentarc, currentparameter,
+                               pointonarc, QuadSurf, OtherQuad, aVertTol);
 	  if (linefound) {
 	    
 #if 1 
@@ -493,14 +497,14 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 	      // deuxieme surface
 	      
 	      goon = PointOnSecondDom (listpnt, Domain, QuadSurf, Normale, 
-				       Vtgint, lin, Done, i);
+                                       Vtgint, lin, Done, i, aVertTol);
 	    }
 	    
 	    if (goon) {
 	      //-- Modification du 4 avril 97    tolerance->Tolarc
 	      //-- on replace sur le vertex la tolerance d entree et 
 	      //-- non la tolerance qui a servi au FindLine
-	      solpnt.SetValue(Psurf,Tolarc,Standard_False);
+              solpnt.SetValue(Psurf, aVertTol, Standard_False);
 	      
 	      U1 = p2d.X(); V1 = p2d.Y();
 	      OtherQuad.Parameters(Psurf,U2,V2);
@@ -513,7 +517,6 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 		Recadre(S1,S2,U2,V2,U1,V1);
 		solpnt.SetParameters(U2,V2,U1,V1);
 	      }
-	      solpnt.SetParameter(paraint);
 	      
 	      if (! currentpointonrst.IsNew()) {
 		vtx = currentpointonrst.Vertex();
@@ -532,12 +535,21 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 	      }
 	      solpnt.SetArc(OnFirst,currentarc, currentparameter,
 			    Transline,Transarc);
-	      if (TheType == IntPatch_Analytic) {
+
-		Handle(IntPatch_ALine)::DownCast (lin)->AddVertex(solpnt);
+              for (TColStd_ListIteratorOfListOfReal anItr(aLParams);
-	      }
+                   anItr.More(); anItr.Next())
-	      else {
+              {
-		Handle(IntPatch_GLine)::DownCast (lin)->AddVertex(solpnt);
+                solpnt.SetParameter(anItr.Value());
-	      }
+                if (TheType == IntPatch_Analytic)
+                {
+                  Handle(IntPatch_ALine)::DownCast(lin)->AddVertex(solpnt);
+                }
+                else
+                {
+                  Handle(IntPatch_GLine)::DownCast(lin)->AddVertex(solpnt);
+                }
+              }
+
 	      Done(i) = 1;
 	      
 	      if (goon) {
@@ -550,6 +562,7 @@ void PutPointsOnLine(const Handle(Adaptor3d_HSurface)& S1,
 		      }		      
 		      else if (Domain->Identical(vtx, vtxbis)) {
 			solpnt.SetVertex(OnFirst,vtxbis);
+                        solpnt.SetTolerance(Tolarc);
 			currentarc = currentpointonrst.Arc();
 			currentparameter = currentpointonrst.Parameter();
 			
@@ -599,7 +612,8 @@ Standard_Boolean  MultiplePoint (const IntPatch_SequenceOfPathPointOfTheSOnBound
                                  TColStd_Array1OfInteger& Done,
                                  TColStd_Array1OfInteger& UsedLine,
                                  const Standard_Integer Index,
-		                 const Standard_Boolean OnFirst) {
+		                 const Standard_Boolean OnFirst,
+                                 const Standard_Real theToler) {
 
 // Traitement des points "multiples".
 
@@ -696,7 +710,7 @@ Standard_Boolean  MultiplePoint (const IntPatch_SequenceOfPathPointOfTheSOnBound
 
 	  intpt.SetArc(OnFirst,currentarc,currentparameter,
 		       Transline,Transarc);
-
+          intpt.SetTolerance(theToler);
 
 	  if (TheType == IntPatch_Analytic) {
 	    Handle(IntPatch_ALine)::DownCast (slinValueii)->Replace(jj,intpt);
@@ -728,6 +742,7 @@ Standard_Boolean  MultiplePoint (const IntPatch_SequenceOfPathPointOfTheSOnBound
 		    }
 	            intpt.SetArc(OnFirst,currentarc,currentparameter,
                                  Transline,Transarc);
+                    intpt.SetTolerance(theToler);
 	            if (TheType == IntPatch_Analytic) {
 	              Handle(IntPatch_ALine)::DownCast (slinValueii)->AddVertex(intpt);
 	            }
@@ -770,7 +785,8 @@ Standard_Boolean PointOnSecondDom (const IntPatch_SequenceOfPathPointOfTheSOnBou
 				   const gp_Vec&    Vtgint,
 				   const Handle(IntPatch_Line)& lin,
 				   TColStd_Array1OfInteger& Done,
-				   const Standard_Integer Index)
+				   const Standard_Integer Index,
+                                   const Standard_Real theToler)
      
 
 // Duplication des points sur domaine de l autre surface.
@@ -841,6 +857,8 @@ Standard_Boolean PointOnSecondDom (const IntPatch_SequenceOfPathPointOfTheSOnBou
 	}
 	intpt.SetArc(Standard_False,currentarc,currentparameter,
 	             Transline,Transarc);
+        intpt.SetTolerance(theToler);
+
 	if (TheType == IntPatch_Analytic) {
 	  Handle(IntPatch_ALine)::DownCast (lin)->Replace(jj,intpt);
 	}
@@ -871,6 +889,7 @@ Standard_Boolean PointOnSecondDom (const IntPatch_SequenceOfPathPointOfTheSOnBou
 		  }
 	          intpt.SetArc(Standard_False,currentarc,currentparameter,
                                Transline,Transarc);
+                  intpt.SetTolerance(theToler);
 	          if (TheType == IntPatch_Analytic) {
 	            Handle(IntPatch_ALine)::DownCast (lin)->AddVertex(intpt);
 	          }
@@ -905,31 +924,36 @@ Standard_Boolean PointOnSecondDom (const IntPatch_SequenceOfPathPointOfTheSOnBou
 
 
 
-Standard_Boolean FindLine (gp_Pnt& Psurf,
+Standard_Boolean FindLine(gp_Pnt& Psurf,
-			   const IntPatch_SequenceOfLine& slin,
+                          const IntPatch_SequenceOfLine& slin,
-			   const Standard_Real Tol,
+                          const Standard_Real Tol,
-			   Standard_Real& Paraint,
+                          TColStd_ListOfReal& theLParams,
-			   gp_Vec& Vtgtint,
+                          gp_Vec& Vtgtint,
-			   Standard_Integer& Range,
+                          Standard_Integer& theLineIdx,
-			   Standard_Integer OnlyThisLine,
+                          Standard_Integer OnlyThisLine,
-			   const Handle(Adaptor2d_HCurve2d)& thearc,
+                          const Handle(Adaptor2d_HCurve2d)& thearc,
-			   Standard_Real& theparameteronarc,
+                          Standard_Real& theparameteronarc,
-			   gp_Pnt& thepointonarc,
+                          gp_Pnt& thepointonarc,
-			   const IntSurf_Quadric& QuadSurf)
+                          const IntSurf_Quadric& QuadSurf1,
+                          const IntSurf_Quadric& QuadSurf2,
+                          Standard_Real& theOutputToler)
 {
+  if ((QuadSurf1.Distance(Psurf) > Tol) || (QuadSurf2.Distance(Psurf) > Tol))
+    return Standard_False;
 
 // Traitement du point de depart ayant pour representation Psurf
 // dans l espace. On recherche la ligne d intersection contenant ce point.
 // On a en sortie la ligne, et le parametre et sa tangente du point sur
 // la ligne d intersection.
-
+  const Standard_Real aSqTol = Tol*Tol;
-  Standard_Real distmin = RealLast();
+  Standard_Real aSqDistMin = RealLast();
-  Standard_Real dist,para;
+  Standard_Real aSqDist, para;
   Standard_Real lower,upper;
   gp_Pnt pt;
   Standard_Integer i;
   IntPatch_IType typarc;
 
+  Standard_Real aParaInt = RealLast();
   Standard_Integer nblin = slin.Length();
   for (i=1; i<=nblin; i++) {
     if(OnlyThisLine) { i=OnlyThisLine; nblin=0; }
@@ -961,11 +985,12 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 	para = ElCLib::Parameter(Handle(IntPatch_GLine)::DownCast (lin)->Line(),Psurf);
 	if (para <= upper && para >= lower) {
 	  pt = ElCLib::Value(para,Handle(IntPatch_GLine)::DownCast (lin)->Line());
-	  dist = Psurf.Distance(pt);
+          aSqDist = Psurf.SquareDistance(pt);
-	  if (dist< distmin) {
+          if ((aSqDist < aSqTol) && (aSqDist < aSqDistMin))
-	    distmin = dist;
+          {
-	    Paraint = para;
+            aSqDistMin = aSqDist;
-	    Range = i;
+            aParaInt = para;
+            theLineIdx = i;
 	  }
 	}
       }
@@ -977,11 +1002,12 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 	    (para + 2.*M_PI <=upper && para + 2.*M_PI >= lower) ||
 	    (para - 2.*M_PI <=upper && para - 2.*M_PI >= lower)) {
 	  pt = ElCLib::Value(para,Handle(IntPatch_GLine)::DownCast (lin)->Circle());
-	  dist = Psurf.Distance(pt);
+          aSqDist = Psurf.SquareDistance(pt);
-	  if (dist< distmin) {
+          if ((aSqDist < aSqTol) && (aSqDist < aSqDistMin))
-	    distmin = dist;
+          {
-	    Paraint = para;
+            aSqDistMin = aSqDist;
-	    Range = i;
+            aParaInt = para;
+            theLineIdx = i;
 	  }
 	}
       }
@@ -993,11 +1019,12 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 	    (para + 2.*M_PI <=upper && para + 2.*M_PI >= lower) ||
 	    (para - 2.*M_PI <=upper && para - 2.*M_PI >= lower)) {
 	  pt = ElCLib::Value(para,Handle(IntPatch_GLine)::DownCast (lin)->Ellipse());
-	  dist = Psurf.Distance(pt);
+          aSqDist = Psurf.SquareDistance(pt);
-	  if (dist< distmin) {
+          if ((aSqDist < aSqTol) && (aSqDist < aSqDistMin))
-	    distmin = dist;
+          {
-	    Paraint = para;
+            aSqDistMin = aSqDist;
-	    Range = i;
+            aParaInt = para;
+            theLineIdx = i;
 	  }
 	}
       }
@@ -1030,24 +1057,28 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 	    Standard_Real parabis = para+0.0000001;
 	    
 	    pt = ElCLib::Value(para,Parab);
-	    dist = Psurf.Distance(pt);
+            aSqDist = Psurf.SquareDistance(pt);
 	    
-	    gp_Pnt ptbis = ElCLib::Value(parabis,Parab);
+	    const gp_Pnt ptbis = ElCLib::Value(parabis,Parab);
-	    Standard_Real distbis = Psurf.Distance(ptbis);
+	    const Standard_Real distbis = Psurf.Distance(ptbis);
-	    
+            const Standard_Real aDist = Sqrt(aSqDist);
-	    Standard_Real ddist = distbis-dist;
+            const Standard_Real ddist = distbis - aDist;
 	    
 	    //--cout<<" para: "<<para<<"    dist:"<<dist<<"   ddist:"<<ddist<<endl;
 	    
-	    if (dist< distmin) {
+            if ((aSqDist < aSqTol) && (aSqDist < aSqDistMin))
-	      distmin = dist;
+            {
-	      Paraint = para;
+              aSqDistMin = aSqDist;
-	      Range = i;
+              aParaInt = para;
+              theLineIdx = i;
 	    }
-	    if(dist<1.0e-9 && dist>-1.0e-9) { amelioration=100; } 
+            if (aSqDist < Precision::SquarePConfusion())
+            {
+              amelioration = 100;
+            }
 	      
 	    if(ddist>1.0e-9 || ddist<-1.0e-9 ) { 
-	      para=para-dist*(parabis-para)/ddist;
+              para = para - aDist*(parabis - para) / ddist;
 	    }
 	    else { 
 	      amelioration=100;
@@ -1065,11 +1096,12 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 	para = ElCLib::Parameter(Handle(IntPatch_GLine)::DownCast (lin)->Hyperbola(),Psurf);
 	if (para <= upper && para >= lower) {
 	  pt = ElCLib::Value(para,Handle(IntPatch_GLine)::DownCast (lin)->Hyperbola());
-	  dist = Psurf.Distance(pt);
+          aSqDist = Psurf.SquareDistance(pt);
-	  if (dist< distmin) {
+          if ((aSqDist < aSqTol) && (aSqDist < aSqDistMin))
-	    distmin = dist;
+          {
-	    Paraint = para;
+            aSqDistMin = aSqDist;
-	    Range = i;
+            aParaInt = para;
+            theLineIdx = i;
 	  }
 	}
       }
@@ -1077,17 +1109,33 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 
     case IntPatch_Analytic :
       {
-	Handle(IntPatch_ALine) alin (Handle(IntPatch_ALine)::DownCast (lin));
+        Handle(IntPatch_ALine) alin(Handle(IntPatch_ALine)::DownCast(lin));
-	Standard_Boolean FindIsOk = alin->FindParameter(Psurf,para);
+        TColStd_ListOfReal aLParams;
-	if (FindIsOk) {
+        alin->FindParameter(Psurf, aLParams);
-	  pt = alin->Value(para);
+        if (!aLParams.IsEmpty())
-	  dist = Psurf.Distance(pt);
+        {
-	  if (dist< distmin) {
+          // All found distances are already in some internal tolerance
-	    distmin = dist;
+          // set in alin->FindParameter(...) method.
-	    Paraint = para;
+
-	    Range = i;
+          aSqDist = RealLast();
-	  }
+          for (TColStd_ListIteratorOfListOfReal anItr(aLParams);
-	}
+               anItr.More(); anItr.Next())
+          {
+            pt = alin->Value(anItr.Value());
+            const Standard_Real aSqD = Psurf.SquareDistance(pt);
+            if (aSqD < aSqDist)
+            {
+              aSqDist = aSqD;
+            }
+          }
+
+          if (aSqDist < aSqDistMin)
+          {
+            aSqDistMin = aSqDist;
+            theLParams = aLParams;
+            theLineIdx = i;
+          }
+        }
 	else { 
 	  //-- le point n a pas ete trouve par bete projection.
 	  //-- on essaie l intersection avec la restriction en 2d
@@ -1096,19 +1144,20 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
 //	  Standard_Real anpara=para;
 //#endif
 	  gp_Pnt CopiePsurf=Psurf;
-	  Standard_Boolean IntersectIsOk=IntersectionWithAnArc(CopiePsurf,alin,para,thearc,theparamonarc,thepointonarc,QuadSurf,lower,upper,dist);
+          Standard_Boolean IntersectIsOk = IntersectionWithAnArc(CopiePsurf, alin, para, 
-
+                                                                 thearc, theparamonarc,
-	  //--printf("\nIntersectionWithAnArc   %d \n Psurf(%g,%g,%g)->(%g,%g,%g)  dist=%g\n para(%g)->(%g)\n paraonarc(%g)->(%g)",
+                                                                 thepointonarc,
-	  //--   ok,Psurf.X(),Psurf.Y(),Psurf.Z(),thepointonarc.X(),thepointonarc.Y(),thepointonarc.Z(),dist,
+                                                                 QuadSurf1,
-	  //--    anpara,para,theparameteronarc,theparamonarc);
+                                                                 lower, upper);
-	  dist = CopiePsurf.Distance(Psurf);
+          aSqDist = CopiePsurf.SquareDistance(Psurf);
 	  if(IntersectIsOk) {
-	    if(dist<Tol) { 
+            if (aSqDist < aSqTol)
+            {
 	      theparameteronarc = theparamonarc;
 	      Psurf = thepointonarc;
-	      distmin = dist;
+              aSqDistMin = aSqDist;
-	      Paraint = para;
+              theLParams.Append(para);
-	      Range = i;
+              theLineIdx = i;
 	    }
 	  }
 	}
@@ -1124,86 +1173,47 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
     }
   }
 
-  if (distmin > Tol) {
+  if (aSqDistMin == RealLast())
     return Standard_False;
-  }
 
-  typarc = slin.Value(Range)->ArcType();
+  theOutputToler = Max(theOutputToler, Sqrt(aSqDistMin));
 
-  // Calcul de la tangente.
+  typarc = slin.Value(theLineIdx)->ArcType();
+
+  // Computation of tangent vector
   switch (typarc) {
   case IntPatch_Lin :
-    Vtgtint = (*((Handle(IntPatch_GLine)*)&slin(Range)))->Line().Direction();
+    theLParams.Append(aParaInt);
+    Vtgtint = (*((Handle(IntPatch_GLine)*)&slin(theLineIdx)))->Line().Direction();
     break;
   case IntPatch_Circle :
-    Vtgtint = ElCLib::DN(Paraint,(*((Handle(IntPatch_GLine)*)&slin(Range)))->Circle(),1);
+    theLParams.Append(aParaInt);
+    Vtgtint = ElCLib::DN(aParaInt, (*((Handle(IntPatch_GLine)*)&slin(theLineIdx)))->Circle(), 1);
     break;
   case IntPatch_Ellipse :
-    Vtgtint = ElCLib::DN(Paraint,(*((Handle(IntPatch_GLine)*)&slin(Range)))->Ellipse(),1);
+    theLParams.Append(aParaInt);
+    Vtgtint = ElCLib::DN(aParaInt, (*((Handle(IntPatch_GLine)*)&slin(theLineIdx)))->Ellipse(), 1);
     break;
   case IntPatch_Parabola :
-    Vtgtint = ElCLib::DN(Paraint,(*((Handle(IntPatch_GLine)*)&slin(Range)))->Parabola(),1);
+    theLParams.Append(aParaInt);
+    Vtgtint = ElCLib::DN(aParaInt, (*((Handle(IntPatch_GLine)*)&slin(theLineIdx)))->Parabola(), 1);
     break;
   case IntPatch_Hyperbola :
-    Vtgtint = ElCLib::DN(Paraint,(*((Handle(IntPatch_GLine)*)&slin(Range)))->Hyperbola(),1);
+    theLParams.Append(aParaInt);
+    Vtgtint = ElCLib::DN(aParaInt, (*((Handle(IntPatch_GLine)*)&slin(theLineIdx)))->Hyperbola(), 1);
     break;
 
   case IntPatch_Analytic:
     {
-      const Handle(IntPatch_ALine)& alin = (*((Handle(IntPatch_ALine)*)&slin(Range)));
+      if (!Handle(IntPatch_ALine)::DownCast(slin(theLineIdx))->D1(theLParams.Last(), pt, Vtgtint))
-      Standard_Boolean abid = alin->D1(Paraint,pt,Vtgtint);
+      {
-      if (!abid) { 
+        //Previously (before the fix #29807) this code tried to process case
-	Standard_Real domaininf,domainsup,paramproche;
+        //when Handle(IntPatch_ALine)::D1(...) method returns FALSE and
-	Standard_Boolean boolbid;
+        //computed Vtgtint input argument value. Currently, any singularities
-	domaininf = alin->FirstParameter(boolbid);
+        //must be processed by high-level algorithms (IntPatch_SpecialPoints class).
-	domainsup = alin->LastParameter(boolbid);
+        //Therefore this code has been deleted as deprecated.
-	if(Paraint>=domaininf && Paraint<=domainsup) { 
+
-	  Standard_Real DeltaParam = 0.001 * (domainsup-domaininf);
+        Vtgtint.SetCoord(0.0, 0.0, 0.0);
-	  if(Paraint-domaininf >= domainsup-Paraint) {
-	    //-- On decale le point vers le parametre le plus eloigne.
-	    DeltaParam = -DeltaParam;
-	  }
-	  Standard_Integer kountbid = 0;
-	  Standard_Boolean bornok = Standard_True;
-	  paramproche = Paraint;
-	  do { 
-	    paramproche+=DeltaParam;
-	    kountbid++;
-	    gp_Pnt ptbid;
-	    if(paramproche>=domaininf && paramproche<=domainsup) {
-	      abid = alin->D1(paramproche,ptbid,Vtgtint);   
-	    }
-	    else { 
-	      bornok = Standard_False; 
-	    } 
-	  }
-	  while(abid==Standard_False && kountbid<5 && bornok);
-	  //-- Attention aux points de tangence (croisement de 4 lignes )
-	  bornok = Standard_True;
-	  kountbid = 0;
-	  gp_Vec OVtgtint(0.0,0.0,0.0);
-	  paramproche = Paraint;
-	  do { 
-	    paramproche-=DeltaParam;
-	    kountbid++;
-	    gp_Pnt ptbid;
-	    if(paramproche>=domaininf && paramproche<=domainsup) {
-	      abid = alin->D1(paramproche,ptbid,OVtgtint);   
-	    }
-	    else { 
-	      bornok = Standard_False; 
-	    } 
-	  }
-	  while(abid==Standard_False && kountbid<5 && bornok);
-	  if(bornok) { 
-	    paramproche = Vtgtint.Dot(OVtgtint);
-	    if(paramproche<=0.0) abid = Standard_False; 
-	  }
-	}
-	if(!abid) { 
-	  //-- cout << "Pb sur Calcul de derivee 111 " << endl;
-	  Vtgtint.SetCoord(0.,0.,0.);
-	}
       }
     }
     break;
@@ -1218,19 +1228,20 @@ Standard_Boolean FindLine (gp_Pnt& Psurf,
   return Standard_True;
 }
 
-
+//=======================================================================
-Standard_Boolean  SingleLine (const gp_Pnt& Psurf,
+//function : SingleLine
-			      const Handle(IntPatch_Line)& lin,
+//purpose  : Traitement du point de depart ayant pour representation Psurf
-			      const Standard_Real Tol,
+//            dans l espace. On le replace sur la ligne d intersection; On a en sortie
-			      Standard_Real& Paraint,
+//            son parametre et sa tangente sur la ligne d intersection.
-			      gp_Vec& Vtgtint) {		 
+//            La fonction renvoie False si le point projete est a une distance
-
+//            superieure a Tol du point a projeter.
-// Traitement du point de depart ayant pour representation Psurf
+//=======================================================================
-// dans l espace. On le replace sur la ligne d intersection; On a en sortie
+Standard_Boolean  SingleLine(const gp_Pnt& Psurf,
-// son parametre et sa tangente sur la ligne d intersection.
+                             const Handle(IntPatch_Line)& lin,
-// La fonction renvoie False si le point projete est a une distance
+                             const Standard_Real Tol,
-// superieure a Tol du point a projeter.
+                             Standard_Real& Paraint,
-
+                             gp_Vec& Vtgtint)
+{
   IntPatch_IType typarc = lin->ArcType();
 
   Standard_Real parproj = 0.;
@@ -1238,7 +1249,6 @@ Standard_Boolean  SingleLine (const gp_Pnt& Psurf,
   gp_Pnt ptproj;
   Standard_Boolean retvalue;
 
-
   switch (typarc) {
   case IntPatch_Lin :
     parproj = ElCLib::Parameter(Handle(IntPatch_GLine)::DownCast (lin)->Line(),Psurf);
@@ -1262,46 +1272,30 @@ Standard_Boolean  SingleLine (const gp_Pnt& Psurf,
     break;
   case IntPatch_Analytic :
     {
-      Handle(IntPatch_ALine) alin (Handle(IntPatch_ALine)::DownCast (lin));
+      Handle(IntPatch_ALine) alin(Handle(IntPatch_ALine)::DownCast(lin));
-      Standard_Boolean ok = alin->FindParameter(Psurf,parproj);
+      TColStd_ListOfReal aLParams;
-      if (ok) {
+      alin->FindParameter(Psurf, aLParams);
-	gp_Pnt ptbid;
+      if (!aLParams.IsEmpty())
-	Standard_Boolean bid = alin->D1(parproj,ptbid,tgint);
+      {
-	if (!bid) { 
+        ptproj = Psurf;
-	  Standard_Real domaininf,domainsup,paramproche;
+        parproj = aLParams.Last();
-	  Standard_Boolean boolbid;
+        gp_Pnt aPtemp;
-	  domaininf = alin->FirstParameter(boolbid);
+        if (!alin->D1(parproj, aPtemp, tgint))
-	  domainsup = alin->LastParameter(boolbid);
+        {
-	  if(parproj>=domaininf && parproj<=domainsup) { 
+          //Previously (before the fix #29807) this code tried to process case
-	    Standard_Real DeltaParam = 0.001 * (domainsup-domaininf);
+          //when Handle(IntPatch_ALine)::D1(...) method returns FALSE and
-	    if(parproj-domaininf >= domainsup-parproj) {
+          //computed Vtgtint input argument value. Currently, any singularities
-	      //-- On decale le point vers le parametre le plus eloigne.
+          //must be processed by high-level algorithms (IntPatch_SpecialPoints class).
-	      DeltaParam = -DeltaParam;
+          //Therefore this code has been deleted as deprecated.
-	    }
+
-	    Standard_Integer kountbid = 0;
+          tgint.SetCoord(0.0, 0.0, 0.0);
-	    paramproche = parproj;
+        }
-	    do { 
-	      paramproche+=DeltaParam;
-	      kountbid++;
-	      bid = alin->D1(paramproche,ptbid,tgint);   
-	    }
-	    while(bid==Standard_False && kountbid<5);
-	    ptproj = Psurf; 
-	  }
-	  if(!bid) { 
-	    //-- cout << "Pb sur Calcul de derivee ALine " << endl;
-	    tgint.SetCoord(0.,0.,0.);
-	    return(Standard_False);
-	  }
-	}
-	else { 
-	  ptproj = Psurf; 
-	}
       }
-      else {
+      else
-	//-- cout << "---- Pb sur ligne analytique dans SingleLine" << endl;
+      {
-	//-- cout << "     Find Parameter"<<endl;
+        //-- cout << "---- Pb sur ligne analytique dans SingleLine" << endl;
-	return Standard_False;
+        //-- cout << "     Find Parameter"<<endl;
+        return Standard_False;
       }
     }
     break;
@@ -1485,6 +1479,7 @@ void ProcessSegments (const IntPatch_SequenceOfSegmentOfTheSOnBounds& listedg,
 	  if (EdgeDegenere==Standard_False && dofirst) {
 	    if (ptvtx.Value().Distance(PStartf.Value()) <=TolArc) {
 	      ptvtx.SetMultiple(Standard_True);
+              ptvtx.SetTolerance(TolArc);
 	      if (typ == IntPatch_Analytic) {
 		Handle(IntPatch_ALine)::DownCast (slinj)->Replace(k,ptvtx);
 	      }
@@ -1549,6 +1544,7 @@ void ProcessSegments (const IntPatch_SequenceOfSegmentOfTheSOnBounds& listedg,
 	  if (EdgeDegenere==Standard_False && dolast) {
 	    if (ptvtx.Value().Distance(PStartl.Value()) <=TolArc) {
 	      ptvtx.SetMultiple(Standard_True);
+              ptvtx.SetTolerance(TolArc);
 	      if (typ == IntPatch_Analytic) {
 		Handle(IntPatch_ALine)::DownCast (slinj)->Replace(k,ptvtx);
 	      }
@@ -1988,6 +1984,7 @@ void ProcessRLine (IntPatch_SequenceOfLine& slin,
 	      }
 	      if (keeppoint) {
 		Ptvtx.SetMultiple(Standard_True);
+                Ptvtx.SetTolerance(_TolArc);
 		newptvtx.SetMultiple(Standard_True);
 		
 		if (typ2 == IntPatch_Analytic) {

src/IntPatch/IntPatch_ImpImpIntersection_2.gxx

@@ -307,11 +307,15 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
       }
       nbpt = solrst.NbPoints();
       nbseg= solrst.NbSegments();
-      for (i=1; i<= nbpt; i++) {
+      for (i = 1; i <= nbpt; i++)
-        pnt1.Append(solrst.Point(i));
+      {
+        const IntPatch_ThePathPointOfTheSOnBounds& aPt = solrst.Point(i);
+        pnt1.Append(aPt);
       }
-      for (i=1; i<= nbseg; i++) {
+      for (i = 1; i <= nbseg; i++)
-        edg1.Append(solrst.Segment(i));
+      {
+        const IntPatch_TheSegmentOfTheSOnBounds& aSegm = solrst.Segment(i);
+        edg1.Append(aSegm);
       }
       nosolonS1 = (nbpt == 0) && (nbseg == 0);
 
@@ -335,15 +339,19 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
       if (solrst.AllArcSolution() && typs1 == typs2) {
         all2 = Standard_True;
       }
+
       nbpt = solrst.NbPoints();
       nbseg= solrst.NbSegments();
       for (i=1; i<= nbpt; i++) {
-        pnt2.Append(solrst.Point(i));
+        const IntPatch_ThePathPointOfTheSOnBounds& aPt = solrst.Point(i);
+        pnt2.Append(aPt);
       }
     
       for (i=1; i<= nbseg; i++) {
-        edg2.Append(solrst.Segment(i));
+        const IntPatch_TheSegmentOfTheSOnBounds& aSegm = solrst.Segment(i);
+        edg2.Append(aSegm);
       }
+      
       nosolonS2 = (nbpt == 0) && (nbseg == 0);
 
       if (nosolonS2 && all2) {  // cas de face sans restrictions

src/IntPatch/IntPatch_ImpImpIntersection_4.gxx

@@ -29,12 +29,10 @@ static void ShortCosForm( const Standard_Real theCosFactor,
                           Standard_Real& theCoeff,
                           Standard_Real& theAngle);
 //
-static 
+static Standard_Boolean ExploreCurve(const gp_Cone& theCo,
-  Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
+                                     IntAna_Curve& aC,
-				const gp_Cone& aCo,
+                                     const Standard_Real aTol,
-				IntAna_Curve& aC,
+                                     IntAna_ListOfCurve& aLC);
-				const Standard_Real aTol,
-				IntAna_ListOfCurve& aLC);
 
 static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
                                       const Standard_Real theUlTarget,
@@ -511,13 +509,6 @@ private:
   const Standard_Boolean myIsReverse;
 };
 
-static 
-  Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
-				const gp_Cone& aCo,
-				IntAna_Curve& aC,
-				const Standard_Real aTol,
-				IntAna_ListOfCurve& aLC);
-
 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const IntSurf_Quadric& theQuad2,
                                   const Handle(IntSurf_LineOn2S)& theLine,
@@ -840,6 +831,7 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
 	if (!procf) {
 	  d=ptf.Distance(ptsol.Value());
 	  if (d <= Tol) {
+            ptsol.SetTolerance(Tol);
 	    if (!ptsol.IsMultiple()) {
 	      //-- le point ptsol (de aligold) est declare multiple sur aligold
 	      Multpoint = Standard_True;
@@ -858,6 +850,7 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
 	}
 	if (!procl) {
 	  if (ptl.Distance(ptsol.Value()) <= Tol) {
+            ptsol.SetTolerance(Tol);
 	    if (!ptsol.IsMultiple()) {
 	      Multpoint = Standard_True;
 	      ptsol.SetMultiple(Standard_True);
@@ -888,6 +881,8 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
       }
     }
   }
+  
+  ptsol.SetTolerance(Tol);
   if (!procf && !procl) {
     Quad1.Parameters(ptf,U1,V1);
     Quad2.Parameters(ptf,U2,V2);
@@ -4136,7 +4131,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
 	//curvsol = anaint.Curve(i);
 	aC=anaint.Curve(i);
 	aLC.Clear();
-	ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
+	ExploreCurve(Co, aC, 10.*Tol, aLC);
 	//
 	aIt.Initialize(aLC);
 	for (; aIt.More(); aIt.Next()) {
@@ -4209,61 +4204,69 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
 }
 //=======================================================================
 //function : ExploreCurve
-//purpose  : 
+//purpose  : Splits aC on several curves in the cone apex points.
 //=======================================================================
-Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
+Standard_Boolean ExploreCurve(const gp_Cone& theCo,
-			      const gp_Cone& aCo,
+                              IntAna_Curve& theCrv,
-			      IntAna_Curve& aC,
+                              const Standard_Real theTol,
-			      const Standard_Real aTol,
+                              IntAna_ListOfCurve& theLC)
-			      IntAna_ListOfCurve& aLC)
-			      
 {
-  Standard_Boolean bFind=Standard_False;
+  const Standard_Real aSqTol = theTol*theTol;
-  Standard_Real aTheta, aT1, aT2, aDst;
+  const gp_Pnt aPapx(theCo.Apex());
-  gp_Pnt aPapx, aPx;
+
+  Standard_Real aT1, aT2;
+  theCrv.Domain(aT1, aT2);
+
+  theLC.Clear();
   //
-  //aC.Dump();
+  TColStd_ListOfReal aLParams;
-  //
+  theCrv.FindParameter(aPapx, aLParams);
-  aLC.Clear();
+  if (aLParams.IsEmpty())
-  aLC.Append(aC);
+  {
-  //
+    theLC.Append(theCrv);
-  aPapx=aCo.Apex();
+    return Standard_False;
-  //
-  aC.Domain(aT1, aT2);
-  //
-  aPx=aC.Value(aT1);
-  aDst=aPx.Distance(aPapx);
-  if (aDst<aTol) {
-    return bFind;
   }
-  aPx=aC.Value(aT2);
+
-  aDst=aPx.Distance(aPapx);
+  for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
-  if (aDst<aTol) {
+  {
-    return bFind;
+    Standard_Real aPrm = anItr.Value();
+
+    if ((aPrm - aT1) < Precision::PConfusion())
+      continue;
+
+    Standard_Boolean isLast = Standard_False;
+    if ((aT2 - aPrm) < Precision::PConfusion())
+    {
+      aPrm = aT2;
+      isLast = Standard_True;
+    }
+
+    const gp_Pnt aP = theCrv.Value(aPrm);
+    const Standard_Real aSqD = aP.SquareDistance(aPapx);
+    if (aSqD < aSqTol)
+    {
+      IntAna_Curve aC1 = theCrv;
+      aC1.SetDomain(aT1, aPrm);
+      aT1 = aPrm;
+      theLC.Append(aC1);
+    }
+
+    if (isLast)
+      break;
   }
-  //
+
-  bFind=aC.FindParameter(aPapx, aTheta);
+  if (theLC.IsEmpty())
-  if (!bFind){
+  {
-    return bFind;
+    theLC.Append(theCrv);
+    return Standard_False;
   }
-  //
+
-  aPx=aC.Value(aTheta);
+  if ((aT2 - aT1) > Precision::PConfusion())
-  aDst=aPx.Distance(aPapx);
+  {
-  if (aDst>aTol) {
+    IntAna_Curve aC1 = theCrv;
-    return !bFind;
+    aC1.SetDomain(aT1, aT2);
+    theLC.Append(aC1);
   }
-  //
+
-  // need to be splitted at aTheta
+  return Standard_True;
-  IntAna_Curve aC1, aC2;
-  //
-  aC1=aC;
-  aC1.SetDomain(aT1, aTheta);
-  aC2=aC;
-  aC2.SetDomain(aTheta, aT2);
-  //
-  aLC.Clear();
-  aLC.Append(aC1);
-  aLC.Append(aC2);
-  //
-  return bFind;
 }

src/IntPatch/IntPatch_ImpPrmIntersection.cxx

@@ -19,6 +19,7 @@
 #include <Adaptor3d_HSurface.hxx>
 #include <Adaptor3d_TopolTool.hxx>
 #include <ElCLib.hxx>
+#include <ElSLib.hxx>
 #include <IntPatch_ArcFunction.hxx>
 #include <IntPatch_PointLine.hxx>
 #include <IntPatch_RLine.hxx>
@@ -87,6 +88,7 @@ static IntPatch_SpecPntType IsSeamOrPole(const Handle(Adaptor3d_HSurface)& theQS
                                          const Handle(IntSurf_LineOn2S)& theLine,
                                          const Standard_Boolean IsReversed,
                                          const Standard_Integer theRefIndex,
+                                         const Standard_Real theTol3D,
                                          const Standard_Real theDeltaMax)
 {
   if((theRefIndex < 1) || (theRefIndex >= theLine->NbPoints()))
@@ -96,6 +98,8 @@ static IntPatch_SpecPntType IsSeamOrPole(const Handle(Adaptor3d_HSurface)& theQS
   Standard_Real aUQRef, aVQRef, aUPRef, aVPRef;
   Standard_Real aUQNext, aVQNext, aUPNext, aVPNext;
 
+  const gp_Pnt &aP3d = theLine->Value(theRefIndex + 1).Value();
+
   if(IsReversed)
   {
     theLine->Value(theRefIndex).Parameters  (aUPRef, aVPRef, aUQRef, aVQRef);
@@ -109,6 +113,28 @@ static IntPatch_SpecPntType IsSeamOrPole(const Handle(Adaptor3d_HSurface)& theQS
 
   const GeomAbs_SurfaceType aType = theQSurf->GetType();
 
+  if ((aType == GeomAbs_Cone) && 
+      (theQSurf->Cone().Apex().SquareDistance(aP3d) < theTol3D*theTol3D))
+  {
+    return IntPatch_SPntPoleSeamU;
+  }
+  else if (aType == GeomAbs_Sphere)
+  {
+    const Standard_Real aSqTol = theTol3D*theTol3D;
+    gp_Pnt aP(ElSLib::Value(0.0, M_PI_2, theQSurf->Sphere()));
+    if (aP.SquareDistance(aP3d) < aSqTol)
+    {
+      return IntPatch_SPntPoleSeamU;
+    }
+    
+    aP = ElSLib::Value(0.0, -M_PI_2, theQSurf->Sphere());
+    if (aP.SquareDistance(aP3d) < aSqTol)
+    {
+      return IntPatch_SPntPoleSeamU;
+    }
+  }
+  
+
   const Standard_Real aDeltaU = Abs(aUQRef - aUQNext);
 
   if((aType != GeomAbs_Torus) && (aDeltaU < theDeltaMax))
@@ -2628,6 +2654,15 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
                                                               PrePoint, IsReversed))
       {
         sline->Add(PrePoint);
+
+        //Avoid adding duplicate points.
+        for (;aFindex <= aLindex; aFindex++)
+        {
+          if (!PrePoint.IsSame(aSSLine->Value(aFindex), theTolTang))
+          {
+            break;
+          }
+        }
       }
       else
       {
@@ -2658,7 +2693,8 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
       DetectOfBoundaryAchievement(theQSurf, IsReversed, aSSLine,
                                   k, aTOL2D, sline, isOnBoundary);
 
-      aPrePointExist = IsSeamOrPole(theQSurf, aSSLine, IsReversed, k - 1, aDeltaUmax);
+      aPrePointExist = IsSeamOrPole(theQSurf, aSSLine, IsReversed,
+                                    k - 1, theTolTang, aDeltaUmax);
 
       if (isOnBoundary && (aPrePointExist != IntPatch_SPntPoleSeamU))
       {
@@ -2742,7 +2778,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
           aSupBound(3) = theQSurf->LastUParameter();
 
           IntPatch_SpecialPoints::
-                      AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_False,
+                      AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_False, 0.0,
                                         aTol, aStartPoint, anInfBound, aSupBound,
                                         aNewPoint, IsReversed);
         }
@@ -2752,9 +2788,10 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
 
           IntPatch_Point aVert;
           aVert.SetValue(aRefPt);
+          aVert.SetTolerance(theTolTang);
 
           if(IntPatch_SpecialPoints::
-                      AddSingularPole(theQSurf, thePSurf, aRefPt, theTolTang,
+                      AddSingularPole(theQSurf, thePSurf, aRefPt,
                                       aVert, aNewPoint, IsReversed))
           {
             aPrePointExist = IntPatch_SPntPole;
@@ -2823,7 +2860,7 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
           aSupBound(3) = theQSurf->LastVParameter();
 
           IntPatch_SpecialPoints::
-                AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_True, aTol,
+                AddPointOnUorVIso(theQSurf, thePSurf, aRefPt, Standard_True, 0.0, aTol,
                                   aStartPoint, anInfBound, aSupBound, aNewPoint,
                                   IsReversed);
         }
@@ -2878,7 +2915,9 @@ static Standard_Boolean DecomposeResult(const Handle(IntPatch_PointLine)& theLin
     if(sline->NbPoints() == 1)
     {
       flNextLine = Standard_True;
-      aFindex = aBindex;
+      
+      if (aFindex < aBindex)
+        aFindex = aBindex;
 
       //Go to the next part of aSSLine
       //because we cannot create the line

src/IntPatch/IntPatch_Intersection.cxx

@@ -1432,8 +1432,32 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the
                                 theS1->IsVPeriodic()? theS1->VPeriod() : 0.0,
                                 theS2->IsUPeriodic()? theS2->UPeriod() : 0.0,
                                 theS2->IsVPeriodic()? theS2->VPeriod() : 0.0};
+
+    NCollection_List<gp_Pnt> aListOfCriticalPoints;
+
+    if (theS1->GetType() == GeomAbs_Cone)
+    {
+      aListOfCriticalPoints.Append(theS1->Cone().Apex());
+    }
+    else if (theS1->GetType() == GeomAbs_Sphere)
+    {
+      aListOfCriticalPoints.Append(theS1->Value(0.0, M_PI_2));
+      aListOfCriticalPoints.Append(theS1->Value(0.0, -M_PI_2));
+    }
+
+    if (theS2->GetType() == GeomAbs_Cone)
+    {
+      aListOfCriticalPoints.Append(theS2->Cone().Apex());
+    }
+    else if (theS2->GetType() == GeomAbs_Sphere)
+    {
+      aListOfCriticalPoints.Append(theS2->Value(0.0, M_PI_2));
+      aListOfCriticalPoints.Append(theS2->Value(0.0, -M_PI_2));
+    }
+
     IntPatch_WLineTool::ExtendTwoWLines(slin, theS1, theS2, TolTang,
-                                        anArrOfPeriod, aBx1, aBx2);
+                                        anArrOfPeriod, aBx1, aBx2,
+                                        aListOfCriticalPoints);
   }
 }
 

src/IntPatch/IntPatch_Point.cxx

@@ -153,18 +153,18 @@ void IntPatch_Point::Dump() const {
   Standard_Real u1,v1,u2,v2;
   pt.Parameters(u1,v1,u2,v2);
 
-  printf("P(%15.10f,%15.10f,%15.10f) UV1(%15.10f,%15.10f)  UV2(%15.10f,%15.10f) (Para:%15.10f)\n",
+  printf("P(%+10.20f,%+10.20f,%+10.20f) UV1(%+10.20f,%+10.20f)  UV2(%+10.20f,%+10.20f) (Para:%+10.20f)\n",
 	 (double)(pt.Value().X()),
 	 (double)(pt.Value().Y()),
 	 (double)(pt.Value().Z()),
 	 (double)u1,(double)v1,(double)u2,(double)v2,(double)para);
   if(onS1)
-    printf("*OnS1*  par=%15.10f arc1=%10p", (double)prm1, (void*)arcS1.operator->());
+    printf("*OnS1*  par=%+10.20f arc1=%10p", (double)prm1, (void*)arcS1.operator->());
   if(vtxonS1)
     printf(" *Vtx1*  vtx1=%10p", (void*)vS1.operator->());
   if(onS1 || vtxonS1) printf("\n");
   if(onS2)
-    printf("*OnS2*  par=%15.10f arc2=%10p", (double)prm2, (void*)arcS2.operator->());
+    printf("*OnS2*  par=%+10.20f arc2=%10p", (double)prm2, (void*)arcS2.operator->());
   if(vtxonS2)
     printf(" *Vtx2*  vtx2=%10p", (void*)vS2.operator->());
 
@@ -184,4 +184,6 @@ void IntPatch_Point::Dump() const {
     }
     cout<<endl;
   }
+  cout << "\nToler = " << Tolerance();
+  cout << endl;
 }

src/IntPatch/IntPatch_Point.hxx

@@ -222,7 +222,10 @@ private:
   Standard_Real para;
   Standard_Real tol;
   Standard_Boolean tgt;
+
+  //TRUE if the point is in several intersection lines
   Standard_Boolean mult;
+
   Standard_Boolean onS1;
   Standard_Boolean vtxonS1;
   Handle(Adaptor3d_HVertex) vS1;

src/IntPatch/IntPatch_PointLine.cxx

@@ -119,8 +119,13 @@ Standard_Real IntPatch_PointLine::
   const gp_Vec aCTan(aN1.Crossed(aN2));
   const Standard_Real aSqMagnFDer = aCTan.SquareMagnitude();
   
-  if(aSqMagnFDer < aSqSmallValue)
+  if (aSqMagnFDer < 1.0e-8)
+  {
+    // Use 1.0e-4 (instead of aSmallValue) to provide
+    // stable computation between different platforms.
+    // See test bugs modalg_7 bug29807_sc01
     return -1.0;
+  }
 
   Standard_Real aDuS1 = 0.0, aDvS1 = 0.0, aDuS2 = 0.0, aDvS2 = 1.0;
 
@@ -160,7 +165,8 @@ Standard_Real IntPatch_PointLine::
   const Standard_Real aDetSyst = aB*aB - aA*aC;
 
   if(Abs(aDetSyst) < aSmallValue)
-  {//Indetermined system solution
+  {
+    //Undetermined system solution
     return -1.0;
   }
 

src/IntPatch/IntPatch_SpecialPoints.cxx

@@ -16,6 +16,7 @@
 #include <IntPatch_SpecialPoints.hxx>
 
 #include <Adaptor3d_HSurface.hxx>
+#include <ElCLib.hxx>
 #include <Extrema_ExtPS.hxx>
 #include <Extrema_GenLocateExtPS.hxx>
 #include <Geom_ConicalSurface.hxx>
@@ -35,10 +36,12 @@ class FuncPreciseSeam: public math_FunctionSetWithDerivatives
 public:
   FuncPreciseSeam(const Handle(Adaptor3d_HSurface)& theQSurf, // quadric
                   const Handle(Adaptor3d_HSurface)& thePSurf, // another surface
-                  const Standard_Boolean isTheUSeam): 
+                  const Standard_Boolean isTheUSeam,
+                  const Standard_Real theIsoParameter): 
         myQSurf(theQSurf),
         myPSurf(thePSurf),
-        mySeamCoordInd(isTheUSeam? 1 : 0) // Defines, U- or V-seam is used
+        mySeamCoordInd(isTheUSeam? 1 : 0), // Defines, U- or V-seam is used
+        myIsoParameter(theIsoParameter)
   {
   };
   
@@ -58,7 +61,7 @@ public:
     try
     {
       const Standard_Integer anIndX = theX.Lower(), anIndF = theF.Lower();
-      Standard_Real aUV[] = {0.0, 0.0};
+      Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
       aUV[mySeamCoordInd] = theX(anIndX+2);
       const gp_Pnt aP1(myPSurf->Value(theX(anIndX), theX(anIndX+1)));
       const gp_Pnt aP2(myQSurf->Value(aUV[0], aUV[1]));
@@ -81,7 +84,7 @@ public:
       const Standard_Integer anIndX = theX.Lower(),
                              anIndRD = theD.LowerRow(),
                              anIndCD = theD.LowerCol();
-      Standard_Real aUV[] = {0.0, 0.0};
+      Standard_Real aUV[] = {myIsoParameter, myIsoParameter};
       aUV[mySeamCoordInd] = theX(anIndX+2);
 
       gp_Pnt aPt;
@@ -131,10 +134,30 @@ private:
   const Handle(Adaptor3d_HSurface)& myQSurf;
   const Handle(Adaptor3d_HSurface)& myPSurf;
 
-  // 1 for U-coordinate, 0 - for V one. 
+  //! 1 for U-coordinate, 0 - for V one. 
   const Standard_Integer mySeamCoordInd;
+
+  //! Constant parameter of iso-line
+  const Standard_Real myIsoParameter;
 };
 
+//=======================================================================
+//function : GetTangent
+//purpose  : Computes tangent having the given parameter.
+//            See calling method(s) for detailed information
+//=======================================================================
+static inline void GetTangent(const Standard_Real theConeSemiAngle,
+                              const Standard_Real theParameter,
+                              gp_XYZ& theResult)
+{
+  const Standard_Real aW2 = theParameter*theParameter;
+  const Standard_Real aCosUn = (1.0 - aW2) / (1.0 + aW2);
+  const Standard_Real aSinUn = 2.0*theParameter / (1.0 + aW2);
+
+  const Standard_Real aTanA = Tan(theConeSemiAngle);
+  theResult.SetCoord(aTanA*aCosUn, aTanA*aSinUn, 1.0);
+}
+
 //=======================================================================
 //function : IsPointOnSurface
 //purpose  : Checks if thePt is in theSurf (with given tolerance).
@@ -281,6 +304,7 @@ Standard_Boolean IntPatch_SpecialPoints::
                                         const Handle(Adaptor3d_HSurface)& thePSurf,
                                         const IntSurf_PntOn2S& theRefPt,
                                         const Standard_Boolean theIsU,
+                                        const Standard_Real theIsoParameter,
                                         const math_Vector& theToler,
                                         const math_Vector& theInitPoint,
                                         const math_Vector& theInfBound,
@@ -292,7 +316,7 @@ Standard_Boolean IntPatch_SpecialPoints::
   IntSurf::SetPeriod(theIsReversed ? thePSurf : theQSurf,
                      theIsReversed ? theQSurf : thePSurf, anArrOfPeriod);
 
-  FuncPreciseSeam aF(theQSurf, thePSurf, theIsU);
+  FuncPreciseSeam aF(theQSurf, thePSurf, theIsU, theIsoParameter);
 
   math_FunctionSetRoot aSRF(aF, theToler);
   aSRF.Perform(aF, theInitPoint, theInfBound, theSupBound);
@@ -323,6 +347,432 @@ Standard_Boolean IntPatch_SpecialPoints::
   return Standard_True;
 }
 
+//=======================================================================
+//function : ProcessSphere
+//purpose  :
+/*
+The intersection point (including the pole)
+must be satisfied to the following system:
+
+    \left\{\begin{matrix}
+    R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
+    R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
+    R*\sin (V_{q})=S_{z}(U_{s},V_{s})
+    \end{matrix}\right,
+where
+  R is the radius of the sphere;
+  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
+  @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
+  @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
+
+Consequently (from first two equations),
+  \left\{\begin{matrix}
+  \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
+  \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
+  \end{matrix}\right.
+
+For pole,
+  V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
+
+Therefore, computation U_{q} directly is impossibly.
+
+Let @V_{q}@ tends to @\pm \pi /2@.
+Then (indeterminate form is evaluated in accordance of L'Hospital rule),
+  \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)}
+  \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}=
+  -\lim_{V_{q} \to (\pi /2-0)}
+  \frac{\frac{\partial S_{x}}
+  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{q}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =
+  -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
+  {\partial U_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} U_{s}}) =
+  -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
+  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}}
+  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
+  {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}).
+
+Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
+
+Let mean, that
+  \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
+  \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
+
+From the 3rd equation of the system, we obtain
+  \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
+  \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
+or
+  R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
+  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
+
+If @V_{q}=\pm \pi /2@, then
+  \frac{\partial S_{z}}{\partial U_{s}}*
+  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
+  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
+
+Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
+  \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
+  -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
+  {\frac{\partial S_{z}}{\partial U_{s}}}.
+
+If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
+  \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
+  -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
+  {\frac{\partial S_{z}}{\partial V_{s}}}
+
+Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} =
+\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
+The reason is written below.
+*/
+//=======================================================================
+Standard_Boolean IntPatch_SpecialPoints::ProcessSphere(const IntSurf_PntOn2S& thePtIso,
+                                                       const gp_Vec& theDUofPSurf,
+                                                       const gp_Vec& theDVofPSurf,
+                                                       const Standard_Boolean theIsReversed,
+                                                       const Standard_Real theVquad,
+                                                       Standard_Real& theUquad,
+                                                       Standard_Boolean& theIsIsoChoosen)
+{
+  theIsIsoChoosen = Standard_False;
+
+  //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
+  //Ask to pay attention to the fact that this vector is always normalized.
+  gp_Vec2d aV1;
+
+  if ((Abs(theDUofPSurf.Z()) < Precision::PConfusion()) &&
+      (Abs(theDVofPSurf.Z()) < Precision::PConfusion()))
+  {
+    //Example of this case is an intersection of a plane with a sphere
+    //when the plane tangents the sphere in some pole (i.e. only one 
+    //intersection point, not line). In this case, U-coordinate of the
+    //sphere is undefined (can be realy anything).
+    //Another reason is that we have tangent zone around the pole
+    //(see bug #26576).
+    //Computation of correct value of theUquad is impossible.
+    //Therefore, (in order to return something) we will consider
+    //the intersection line goes along some isoline in neighborhood
+    //of the pole.
+
+#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
+    cout << "Cannot find UV-coordinate for quadric in the pole."
+      " See considered comment above. IntPatch_SpecialPoints.cxx,"
+      " ProcessSphere(...)" << endl;
+#endif
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    theUquad = aUIso;
+    theIsIsoChoosen = Standard_True;
+  }
+  else
+  {
+    if (Abs(theDUofPSurf.Z()) > Abs(theDVofPSurf.Z()))
+    {
+      const Standard_Real aDusDvs = theDVofPSurf.Z() / theDUofPSurf.Z();
+      aV1.SetCoord(theDUofPSurf.X()*aDusDvs - theDVofPSurf.X(),
+                   theDUofPSurf.Y()*aDusDvs - theDVofPSurf.Y());
+    }
+    else
+    {
+      const Standard_Real aDvsDus = theDUofPSurf.Z() / theDVofPSurf.Z();
+      aV1.SetCoord(theDVofPSurf.X()*aDvsDus - theDUofPSurf.X(),
+                   theDVofPSurf.Y()*aDvsDus - theDUofPSurf.Y());
+    }
+
+    aV1.Normalize();
+
+    if (Abs(aV1.X()) > Abs(aV1.Y()))
+      theUquad = Sign(asin(aV1.Y()), theVquad);
+    else
+      theUquad = Sign(acos(aV1.X()), theVquad);
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : ProcessCone
+//purpose  : 
+/*
+The intersection point (including the pole)
+must be satisfied to the following system:
+
+  \left\{\begin {matrix}
+  (V_{q}\sin(a) + R)*\cos(U_{q})) = S_{x}(U_{s}, V_{s})\\
+  (V_{q}\sin(a) + R)*\sin(U_{q})) = S_{y}(U_{s}, V_{s})\\
+  V_{q}\cos(a) = S_{z}(U_{s}, V_{s})
+  \end {matrix}\right,
+where 
+  R is the radius of the cone;
+  a is its semi-angle;
+  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
+  @U_{s}@ and @V_{s}@ are parameters on the parametric surface;
+  @U_{q}@ and @V_{q}@ are equal to theUquad and theVquad correspondingly.
+
+Consequently (from first two equations), 
+  \left\{\begin{matrix}
+  \cos(U_{q})=\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}\\
+  \sin(U_{q})=\frac{S_{y}(U_{s}, V_{s})}{(V_{q}\sin(a)+R)}
+  \end{matrix}\right.
+
+For pole, the denominator of these two equations is equal to 0.
+Therefore, computation U_{q} directly is impossibly.
+
+Let @V_{q}@ tends to @\frac{-R}{\sin(a)})@.
+Then (indeterminate form is evaluated in accordance of L'Hospital rule),
+
+  \cos (U_{q}) = 
+  \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{S_{x}(U_{s},V_{s})}{(V_{q}\sin(a)+R)}=
+  \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dU_{s}}{dV_{q}}*
+  (\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}{\partial V_{s}}*
+  \frac{dV_{s}}{dU_{s}})=
+  \frac{1}{\sin(a)}* \lim_{V_{q} \to (\frac{-R}{\sin(a)})}\frac{dV_{s}}{dV_{q}}*
+  (\frac{\partial S_{x}}{\partial U_{s}}*
+  \frac{dU_{s}}{dV_{s}}+\frac{\partial S_{x}}{\partial V_{s}})
+
+Analogically for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
+
+After differentiating 3rd equation of the system, we will obtain
+  \cos(a)=\frac{dS_{z}}{dV_{q}}=\frac{dU_{s}}{dV_{q}}*
+  (\frac{\partial S_{z}}{\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*
+  \frac{dV_{s}}{dU_{s}})
+or
+  \frac{dU_{s}}{dV_{q}}=\frac{\cos(a)}{\frac{\partial S_{z}}{\partial U_{s}}+
+  \frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+After substituting we will obtain
+  \cos (U_{q}) = 
+  \cot(a)*\frac{\frac{\partial S_{x}}{\partial U_{s}}+\frac{\partial S_{x}}
+  {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
+  {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+  \sin (U_{q}) =
+  \cot(a)*\frac{\frac{\partial S_{y}}{\partial U_{s}}+\frac{\partial S_{y}}
+  {\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}{\frac{\partial S_{z}}
+  {\partial U_{s}}+\frac{\partial S_{z}}{\partial V_{s}}*\frac{dV_{s}}{dU_{s}}}
+
+So, we have obtained vector with coordinates {@ \cos (U_{q}) @, @ \sin (U_{q}) @}.
+Ask to pay attention to the fact that this vector is always normalized.
+And after normalization this vector will have coordinates
+  {\cos (U_{q}), \sin (U_{q})} = {dS_{x}, dS_{y}}.Normalized().
+
+It means that we have to compute a tangent to the intersection curve in
+the cone apex point. After that, just take its X- and Y-coordinates.
+
+However, we have to compute derivative @\frac{dV_{s}}{dU_{s}}@ in order
+to compute this vector. In order to find this derivative, we use the
+information about direction of tangent to the intersection curve.
+This tangent will be directed along the cone generatrix obtained by intersection
+of the cone with a plane tangent to 2nd (intersected) surface.
+*/
+//=======================================================================
+Standard_Boolean IntPatch_SpecialPoints::ProcessCone(const IntSurf_PntOn2S& thePtIso,
+                                                     const gp_Vec& theDUofPSurf,
+                                                     const gp_Vec& theDVofPSurf,
+                                                     const gp_Cone& theCone,
+                                                     const Standard_Boolean theIsReversed,
+                                                     Standard_Real& theUquad,
+                                                     Standard_Boolean& theIsIsoChoosen)
+{
+  theIsIsoChoosen = Standard_False;
+
+  // A plane tangent to 2nd (intersected) surface.
+  // Its normal.
+  const gp_XYZ aTgPlaneZ(theDUofPSurf.Crossed(theDVofPSurf).XYZ());
+  const Standard_Real aSqModTg = aTgPlaneZ.SquareModulus();
+  if (aSqModTg < Precision::SquareConfusion())
+  {
+    theIsIsoChoosen = Standard_True;
+  }
+
+  gp_XYZ aTgILine[2];
+  const Standard_Integer aNbTangent = !theIsIsoChoosen? 
+                          GetTangentToIntLineForCone(theCone.SemiAngle(),
+                                                     aTgPlaneZ.Divided(Sqrt(aSqModTg)),
+                                                     aTgILine) : 0;
+
+  if (aNbTangent == 0)
+  {
+    theIsIsoChoosen = Standard_True;
+  }
+  else
+  {
+    const Standard_Real aPeriod = M_PI + M_PI;
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    aUIso = ElCLib::InPeriod(aUIso, 0.0, aPeriod);
+
+    // Sought U-parameter in the apex point
+
+    // For 2 possible parameters value,
+    // one will be chosen which is nearer
+    // to aUIso. Following variables will help to chose.
+    Standard_Real aMinDelta = RealLast();
+    for (Standard_Integer anIdx = 0; anIdx < aNbTangent; anIdx++)
+    {
+      // Vector {@\cos(a), \sin(a)@}
+      gp_Vec2d aVecCS(aTgILine[anIdx].X(), aTgILine[anIdx].Y());
+      const Standard_Real aSqMod = aVecCS.SquareMagnitude();
+      if (aSqMod < Precision::SquareConfusion())
+      {
+        theIsIsoChoosen = Standard_True;
+        break;
+      }
+
+      // Normalize
+      aVecCS.Divide(Sqrt(aSqMod));
+
+      // Angle in range [0, PI/2]
+      Standard_Real anUq = (Abs(aVecCS.X()) < Abs(aVecCS.Y())) ? ACos(Abs(aVecCS.X())) : ASin(Abs(aVecCS.Y()));
+
+      // Convert angles to the range [0, 2*PI]
+      if (aVecCS.Y() < 0.0)
+      {
+        if (aVecCS.X() > 0.0)
+        {
+          anUq = -anUq;
+        }
+        else
+        {
+          anUq += M_PI;
+        }
+      }
+      else if (aVecCS.X() < 0.0)
+      {
+        anUq = M_PI - anUq;
+      }
+
+      //Select the parameter the nearest to aUIso
+      anUq = ElCLib::InPeriod(anUq, 0.0, aPeriod);
+      Standard_Real aDelta = Abs(anUq - aUIso);
+      if (aDelta > M_PI)
+        aDelta = aPeriod - aDelta;
+
+      if (aDelta < aMinDelta)
+      {
+        aMinDelta = aDelta;
+        theUquad = anUq;
+      }
+    }
+  }
+
+  if (theIsIsoChoosen)
+  {
+#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
+    cout << "Cannot find UV-coordinate for quadric in the pole."
+      " IntPatch_AddSpecialPoints.cxx, ProcessCone(...)" << endl;
+#endif
+    theIsIsoChoosen = Standard_True;
+
+    Standard_Real aUIso = 0.0, aVIso = 0.0;
+    if (theIsReversed)
+      thePtIso.ParametersOnS2(aUIso, aVIso);
+    else
+      thePtIso.ParametersOnS1(aUIso, aVIso);
+
+    theUquad = aUIso;
+    return Standard_True;
+  }
+  else
+  {
+    return Standard_True;
+  }
+
+  //return Standard_False;
+}
+
+//=======================================================================
+//function : GetTangentToIntLineForCone
+//purpose  : The following conditions must be satisfied:
+//1. The cone is represented in its canonical form.
+//2. The plane goes through the cone apex and has the normal vector thePlnNormal.
+//3. Vector thePlnNormal has already been normalized
+/*
+Let us enter the new coordinate system where the origin will be in the cone apex
+and axes are the same as in World-Coordinate-System (WCS).
+There the horizontal plane (which is parallel XY-plane) with the origin
+(0, 0, 1) will intersect the cone by the circle with center (0, 0, 1),
+direction {0, 0, 1} and radius tg(a) where a is the cone semi-angle.
+Its equation will be
+\left\{\begin{matrix}
+x(U_{n}) = \tan(a)*\cos(U_{n}) = \tan(a)*\frac{1-\tan^{2}(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
+y(U_{n}) = \tan(a)*\sin (U_{n}) = \tan(a)*\frac{2*\tan(U_{n}/2)}{1+\tan^{2}(U_{n}/2)}\\
+z(U_{n}) = 1
+\end{matrix}\right.
+
+The given plane has (in this coordinate system) location (0, 0, 0) and
+the same normal thePlnNormal=={nx,ny,nz}. Its equation is:
+nx*x+ny*y+nz*z==0
+
+After substitution circle's equation to the plane's equation
+we will obtain a quadratic equation
+aA*w^2 + 2*aB*w + aC = 0.
+*/
+//=======================================================================
+Standard_Integer IntPatch_SpecialPoints::GetTangentToIntLineForCone(const Standard_Real theConeSemiAngle,
+                                                                    const gp_XYZ& thePlnNormal,
+                                                                    gp_XYZ theResult[2])
+{
+  const Standard_Real aNullTol = Epsilon(1.0);
+  const Standard_Real aTanA = Tan(theConeSemiAngle);
+  const Standard_Real aA = thePlnNormal.Z() / aTanA - thePlnNormal.X();
+  const Standard_Real aB = thePlnNormal.Y();
+  const Standard_Real aC = thePlnNormal.Z() / aTanA + thePlnNormal.X();
+
+  if (Abs(aA) < aNullTol)
+  {
+    if (Abs(aB) > aNullTol)
+    {
+      //The plane goes along the cone generatrix.
+      GetTangent(theConeSemiAngle, -aC / (aB + aB), theResult[0]);
+      return 1;
+    }
+
+    //The cone and the plane have only one common point.
+    //It is the cone apex.
+    return 0;
+  }
+
+  //Discriminant of this equation is equal to 
+  Standard_Real aDiscr = thePlnNormal.Z() / Sin(theConeSemiAngle);
+  aDiscr = 1.0 - aDiscr*aDiscr;
+
+  if (Abs(aDiscr) < aNullTol)
+  {
+    //The plane goes along the cone generatrix.
+    // Attention! Mathematically, this cond. is equivalent to 
+    // above processed one (Abs(aA) < aNullTol && (Abs(aB) > aNullTol)).
+    // However, we separate this branch in order to eliminate numerical
+    // instability.
+
+    GetTangent(theConeSemiAngle, -aB / aA, theResult[0]);
+    return 1;
+  }
+  else if (aDiscr > 0.0)
+  {
+    const Standard_Real aRD = Sqrt(aDiscr);
+    GetTangent(theConeSemiAngle, (-aB+aRD)/aA, theResult[0]);
+    GetTangent(theConeSemiAngle, (-aB-aRD)/aA, theResult[1]);
+    return 2;
+  }
+
+  // We will never come here.
+  return 0;
+}
+
 //=======================================================================
 //function : AddSingularPole
 //purpose  : theQSurf is the surface possibly containing special point, 
@@ -333,22 +783,11 @@ Standard_Boolean IntPatch_SpecialPoints::
                       AddSingularPole(const Handle(Adaptor3d_HSurface)& theQSurf,
                                       const Handle(Adaptor3d_HSurface)& thePSurf,
                                       const IntSurf_PntOn2S& thePtIso,
-                                      const Standard_Real theTol,
                                       IntPatch_Point& theVertex,
                                       IntSurf_PntOn2S& theAddedPoint,
                                       const Standard_Boolean theIsReversed,
                                       const Standard_Boolean theIsReqRefCheck)
 {
-  const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0;
-  const Standard_Real aUqPeriod = theQSurf->IsUPeriodic() ? theQSurf->UPeriod() : 0.0;
-  const Standard_Real aVpPeriod = thePSurf->IsVPeriodic() ? thePSurf->VPeriod() : 0.0;
-  const Standard_Real aVqPeriod = theQSurf->IsVPeriodic() ? theQSurf->VPeriod() : 0.0;
-
-  const Standard_Real anArrOfPeriod[4] = {theIsReversed? aUpPeriod : aUqPeriod,
-                                          theIsReversed? aVpPeriod : aVqPeriod,
-                                          theIsReversed? aUqPeriod : aUpPeriod,
-                                          theIsReversed? aVqPeriod : aVpPeriod};
-
   //On parametric
   Standard_Real aU0 = 0.0, aV0 = 0.0;
   //aPQuad is Pole
@@ -379,13 +818,13 @@ Standard_Boolean IntPatch_SpecialPoints::
   }
 
   theQSurf->D0(aUquad, aVquad, aPQuad);
-
+  const Standard_Real aTol = theVertex.Tolerance();
-  if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= theTol*theTol))
+  if (theIsReqRefCheck && (aPQuad.SquareDistance(theVertex.Value()) >= aTol*aTol))
   {
     return Standard_False;
   }
 
-  if(!IsPointOnSurface(thePSurf, aPQuad, theTol, aP0, aU0, aV0))
+  if (!IsPointOnSurface(thePSurf, aPQuad, aTol, aP0, aU0, aV0))
   {
     return Standard_False;
   }
@@ -398,25 +837,20 @@ Standard_Boolean IntPatch_SpecialPoints::
   else
     theAddedPoint.SetValue(0.5*(aP0.XYZ() + aPQuad.XYZ()), aUquad, aVquad, aU0, aV0);
 
-  Standard_Boolean isSame = Standard_False;
+  const Standard_Boolean isSame = theAddedPoint.IsSame(theVertex.PntOn2S(),
-
+                                                       Precision::Confusion());
-  if (theAddedPoint.IsSame(theVertex.PntOn2S(), Precision::Confusion()))
-  {
-    isSame = Standard_True;
-  }
 
   //Found pole does not exist in the Walking-line
   //It must be added there (with correct 2D-parameters)
 
-  //2D-parameters of theparametric surface have already been found (aU0, aV0).
+  //2D-parameters of thePSurf surface have already been found (aU0, aV0).
   //Let find 2D-parameters on the quadric.
 
-  //The algorithm depends on the type of the quadric. Here we consider a Sphere only.
+  //The algorithm depends on the type of the quadric.
-  //Analogical result can be made for another types (e.g. cone, but formulas will
+  //Here we consider a Sphere and cone only.
-  //be different) in case of need.
 
-  //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere
+  //First of all, we need in adjusting thePSurf in the coordinate system of the Sphere/Cone
-  //(in order to make the equation of the sphere maximal simple). However, as it will be
+  //(in order to make its equation maximal simple). However, as it will be
   //shown later, thePSurf is used in algorithm in order to get its derivatives.
   //Therefore, for improving performance, transformation of these vectors is enough
   //(there is no point in transformation of full surface).
@@ -439,146 +873,19 @@ Standard_Boolean IntPatch_SpecialPoints::
 
   if(theQSurf->GetType() == GeomAbs_Sphere)
   {
-    //The intersection point (including the pole)
+    if (!ProcessSphere(thePtIso, aVecDu, aVecDv, theIsReversed, 
-    //must be satisfied to the following system:
+                        aVquad, aUquad, isIsoChoosen))
-    
-    //    \left\{\begin{matrix}
-    //    R*\cos (U_{q})*\cos (V_{q})=S_{x}(U_{s},V_{s})
-    //    R*\sin (U_{q})*\cos (V_{q})=S_{y}(U_{s},V_{s})
-    //    R*\sin (V_{q})=S_{z}(U_{s},V_{s})
-    //    \end{matrix}\right,
-    //where 
-    //  R is the radius of the sphere;
-    //  @S_{x}@, @S_{y}@ and @S_{z}@ are X, Y and Z-coordinates of thePSurf;
-    //  @U_{s}@ and @V_{s}@ are equal to aU0 and aV0 corespondingly;
-    //  @U_{q}@ and @V_{q}@ are equal to aUquad and aVquad corespondingly.
-
-    //Consequently (from first two equations), 
-    //  \left\{\begin{matrix}
-    //  \cos (U_{q}) = \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}
-    //  \sin (U_{q}) = \frac{S_{y}(U_{s},V_{s})}{R*\cos (V_{q})}
-    //  \end{matrix}\right.
-
-    //For pole, 
-    //  V_{q}=\pm \pi /2 \Rightarrow \cos (V_{q}) = 0 (denominator is equal to 0).
-
-    //Therefore, computation U_{q} directly is impossibly.
-    //
-    //Let @V_{q}@ tends to @\pm \pi /2@.
-    //Then (indeterminate form is evaluated in accordance of L'Hospital rule),
-    //  \cos (U_{q}) = \lim_{V_{q} \to (\pi /2-0)} 
-    //  \frac{S_{x}(U_{s},V_{s})}{R*\cos (V_{q})}= 
-    //  -\lim_{V_{q} \to (\pi /2-0)}
-    //  \frac{\frac{\partial S_{x}}
-    //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}} 
-    //  {\mathrm{d} V_{q}}+\frac{\partial S_{x}} 
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}} 
-    //  {\mathrm{d} V_{q}}}{R*\sin (V_{q})} =  
-    //  -\frac{1}{R}*\frac{\mathrm{d} U_{s}}
-    //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
-    //  {\partial U_{s}}+\frac{\partial S_{x}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}
-    //  {\mathrm{d} U_{s}}) =
-    //  -\frac{1}{R}*\frac{\mathrm{d} V_{s}}
-    //  {\mathrm{d} V_{q}}*(\frac{\partial S_{x}} 
-    //  {\partial U_{s}}*\frac{\mathrm{d} U_{s}}
-    //  {\mathrm{d} V_{s}}+\frac{\partial S_{x}}
-    //  {\partial V_{s}}).
-
-    //Analogicaly for @\sin (U_{q})@ (@S_{x}@ is substituted to @S_{y}@).
-
-    //Let mean, that
-    //  \cos (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \cos (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
-    //  \sin (U_{q}) \left | _{V_{q} \to (-\pi /2+0)} = \sin (U_{q}) \left | _{V_{q} \to (\pi /2-0)}
-
-    //From the 3rd equation of the system, we obtain
-    //  \frac{\mathrm{d} (R*\sin (V_{q}))}{\mathrm{d} V_{q}} =
-    //  \frac{\mathrm{d} S_{z}(U_{s},V_{s})}{\mathrm{d} V_{q}}
-    //or
-    //  R*\cos (V_{q}) = \frac{\partial S_{z}}{\partial U_{s}}*
-    //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}}.
-
-    //If @V_{q}=\pm \pi /2@, then
-    //  \frac{\partial S_{z}}{\partial U_{s}}*
-    //  \frac{\mathrm{d} U_{s}} {\mathrm{d} V_{q}}+\frac{\partial S_{z}}
-    //  {\partial V_{s}}*\frac{\mathrm{d} V_{s}}{\mathrm{d} V_{q}} = 0.
-
-    //Consequently, if @\frac{\partial S_{z}}{\partial U_{s}} \neq 0 @ then
-    //  \frac{\mathrm{d} U_{s}}{\mathrm{d} V_{s}} =
-    //  -\frac{\frac{\partial S_{z}}{\partial V_{s}}}
-    //  {\frac{\partial S_{z}}{\partial U_{s}}}.
-
-    //If @ \frac{\partial S_{z}}{\partial V_{s}} \neq 0 @ then
-    //  \frac{\mathrm{d} V_{s}}{\mathrm{d} U_{s}} =
-    //  -\frac{\frac{\partial S_{z}}{\partial U_{s}}}
-    //  {\frac{\partial S_{z}}{\partial V_{s}}}
-
-    //Cases, when @ \frac{\partial S_{z}}{\partial U_{s}} = 
-    //\frac{\partial S_{z}}{\partial V_{s}} = 0 @ are not consider here.
-    //The reason is written below.
-
-    //Vector with {@ \cos (U_{q}) @, @ \sin (U_{q}) @} coordinates.
-    //Ask to pay attention to the fact that this vector is always normalyzed.
-    gp_Vec2d aV1;
-
-    if( (Abs(aVecDu.Z()) < Precision::PConfusion()) &&
-        (Abs(aVecDv.Z()) < Precision::PConfusion()))
     {
-      //Example of this case is an intersection of a plane with a sphere
+      return Standard_False;
-      //when the plane tangents the sphere in some pole (i.e. only one 
-      //intersection point, not line). In this case, U-coordinate of the
-      //sphere is undefined (can be realy anything).
-      //Another reason is that we have tangent zone around the pole
-      //(see bug #26576).
-      //Computation of correct value of aUquad is impossible.
-      //Therefore, (in oreder to return something) we will consider
-      //the intersection line goes along some isoline in neighbourhood
-      //of the pole.
-
-#ifdef INTPATCH_ADDSPECIALPOINTS_DEBUG
-      cout << "Cannot find UV-coordinate for quadric in the pole."
-        " See considered comment above. IntPatch_AddSpecialPoints.cxx,"
-        " AddSingularPole(...)" << endl;
-#endif
-      Standard_Real aUIso = 0.0, aVIso = 0.0;
-      if(theIsReversed)
-        thePtIso.ParametersOnS2(aUIso, aVIso);
-      else
-        thePtIso.ParametersOnS1(aUIso, aVIso);
-
-      aUquad = aUIso;
-      isIsoChoosen = Standard_True;
-    }
-    else
-    {
-      if(Abs(aVecDu.Z()) > Abs(aVecDv.Z()))
-      {
-        const Standard_Real aDusDvs = aVecDv.Z()/aVecDu.Z();
-        aV1.SetCoord( aVecDu.X()*aDusDvs - aVecDv.X(),
-                      aVecDu.Y()*aDusDvs - aVecDv.Y());
-      }
-      else
-      {
-        const Standard_Real aDvsDus = aVecDu.Z()/aVecDv.Z();
-        aV1.SetCoord( aVecDv.X()*aDvsDus - aVecDu.X(),
-                      aVecDv.Y()*aDvsDus - aVecDu.Y());
-      }
-
-      aV1.Normalize();
-
-      if(Abs(aV1.X()) > Abs(aV1.Y()))
-        aUquad = Sign(asin(aV1.Y()), aVquad);
-      else
-        aUquad = Sign(acos(aV1.X()), aVquad);
     }
   }
   else //if(theQSurf->GetType() == GeomAbs_Cone)
   {
-    // This case is not processed. However,
+    if (!ProcessCone(thePtIso, aVecDu, aVecDv, theQSurf->Cone(),
-    // it can be done using the same algorithm
+                      theIsReversed, aUquad, isIsoChoosen))
-    // as for sphere (formulas will be different).
+    {
-    return Standard_False;
+      return Standard_False;
+    }
   }
 
   if(theIsReversed)
@@ -594,6 +901,16 @@ Standard_Boolean IntPatch_SpecialPoints::
 
   if (!isIsoChoosen)
   {
+    Standard_Real anArrOfPeriod[4];
+    if (theIsReversed)
+    {
+      IntSurf::SetPeriod(thePSurf, theQSurf, anArrOfPeriod);
+    }
+    else
+    {
+      IntSurf::SetPeriod(theQSurf, thePSurf, anArrOfPeriod);
+    }
+
     AdjustPointAndVertex(theVertex.PntOn2S(), anArrOfPeriod, theAddedPoint);
   }
   else
@@ -606,7 +923,32 @@ Standard_Boolean IntPatch_SpecialPoints::
 
 //=======================================================================
 //function : ContinueAfterSpecialPoint
-//purpose  : 
+//purpose  : If the last point of the line is the pole of the quadric then
+//            the Walking-line has been broken in this point.
+//           However, new line must start from this point. Here we must
+//            find 2D-coordinates of "this new" point.
+/*
+The inters. line in the neighborhood of the Apex/Pole(s) can be
+approximated by the intersection result of the Cone/Sphere with
+the plane going through the Apex/Pole and being tangent to the
+2nd intersected surface. This intersection result is well known.
+
+In case of sphere, the inters. result is a circle.
+If we go along this circle and across the Pole then U-parameter of
+the sphere (@U_{q}@) will change to +/-PI.
+
+In case of cone, the inters. result is two intersected lines (which
+can be merged to one in a special case when the plane goes along
+some generatrix of the cone). The direction of these lines
+are computed by GetTangentToIntLineForCone(...) method).
+
+When the real (not lines) inters. curve goes through the cone apex then
+two variants are possible:
+a) The tangent line to the inters. curve will be left. In this case
+U-parameter of the cone (@U_{q}@) will be change to +/-PI.
+b) Another line (as inters. result of cone + plane) will tangent
+to the inters. curve. In this case @U_{q}@ must be recomputed.
+*/
 //=======================================================================
 Standard_Boolean IntPatch_SpecialPoints::
                   ContinueAfterSpecialPoint(const Handle(Adaptor3d_HSurface)& theQSurf,
@@ -620,36 +962,63 @@ Standard_Boolean IntPatch_SpecialPoints::
   if(theSPType == IntPatch_SPntNone)
     return Standard_False;
 
-  //If the last point of the line is the pole of the quadric.
-  //In this case, Walking-line has been broken in this point.
-  //However, new line must start from this point. Here we must
-  //find its 2D-coordinates.
-
-  //For sphere and cone, some intersection point is satisfied to the system
-  //  \cos(U_{q}) = S_{x}(U_{s},V_{s})/F(V_{q}) 
-  //  \sin(U_{q}) = S_{y}(U_{s},V_{s})/F(V_{q}) 
-
-  //where 
-  //  @S_{x}@, @S_{y}@ are X and Y-coordinates of thePSurf;
-  //  @U_{s}@ and @V_{s}@ are UV-parameters on thePSurf;
-  //  @U_{q}@ and @V_{q}@ are UV-parameters on theQSurf;
-  //  @F(V_{q}) @ is some function, which value independs on @U_{q}@
-  //              (form of this function depends on the type of the quadric).
-
-  //When we go through the pole/apex, the function @F(V_{q}) @ changes sign.
-  //Therefore, some cases are possible, when only @\cos(U_{q}) @ or
-  //only @ \sin(U_{q}) @ change sign.
-
-  //Consequently, when the line goes throug the pole, @U_{q}@ can be
-  //changed on @\pi /2 @ (but not less).
-    
   if(theNewPoint.IsSame(theRefPt, Precision::Confusion(), theTol2D))
   {
     return Standard_False;
   }
 
-  //Here, in case of pole/apex adding, we forbid "jumping" between two neighbor
+  if ((theSPType == IntPatch_SPntPole) && (theQSurf->GetType() == GeomAbs_Cone))
-  //Walking-point with step greater than pi/4
+  {
+    //Check if the condition b) is satisfied.
+    //Repeat the same steps as in 
+    //IntPatch_SpecialPoints::AddSingularPole(...) method.
+
+    //On parametric
+    Standard_Real aU0 = 0.0, aV0 = 0.0;
+    //On quadric
+    Standard_Real aUquad = 0.0, aVquad = 0.0;
+
+    if (theIsReversed)
+      theNewPoint.Parameters(aU0, aV0, aUquad, aVquad);
+    else
+      theNewPoint.Parameters(aUquad, aVquad, aU0, aV0);
+
+    gp_Pnt aPtemp;
+    gp_Vec aVecDu, aVecDv;
+    thePSurf->D1(aU0, aV0, aPtemp, aVecDu, aVecDv);
+
+    //Transforms parametric surface in coordinate-system of the quadric
+    gp_Trsf aTr;
+    aTr.SetTransformation(theQSurf->Cone().Position());
+
+    //Derivatives of transformed thePSurf
+    aVecDu.Transform(aTr);
+    aVecDv.Transform(aTr);
+
+    Standard_Boolean isIsoChoosen = Standard_False;
+    ProcessCone(theRefPt, aVecDu, aVecDv, theQSurf->Cone(),
+                theIsReversed, aUquad, isIsoChoosen);
+
+    theNewPoint.SetValue(!theIsReversed, aUquad, aVquad);
+  }
+
+  //As it has already been said, in case of going through the Pole/Apex,
+  //U-parameter of the quadric surface will change to +/-PI. This rule has some
+  //exceptions:
+  //1. When 2nd surface has C0-continuity in the point common with the Apex/Pole.
+  //    In this case, the tangent line to the intersection curve after the Apex/Pole
+  //    must be totally recomputed according to the new derivatives of the 2nd surface.
+  //    Currently, it is not implemented but will be able to be done after the
+  //    corresponding demand.
+  //2. The inters. curve has C1 continuity but huge curvature in the point common with 
+  //    the Apex/Pole. Existing inters. algorithm does not allow putting many points
+  //    near to the Apex/Pole in order to cover this "sharp" piece of the inters. curve.
+  //    Therefore, we use adjusting U-parameter of the quadric surface with
+  //    period PI/2 instead of 2PI. It does not have any mathematical idea
+  //    but allows creating WLine with more or less uniform distributed points.
+  //    In other words, we forbid "jumping" between two neighbor Walking-points
+  //    with step greater than PI/4.
+
   const Standard_Real aPeriod = (theSPType == IntPatch_SPntPole)? M_PI_2 : 2.0*M_PI;
 
   const Standard_Real aUpPeriod = thePSurf->IsUPeriodic() ? thePSurf->UPeriod() : 0.0;

src/IntPatch/IntPatch_SpecialPoints.hxx

@@ -24,6 +24,9 @@
 #include <Standard_Handle.hxx>
 
 class Adaptor3d_HSurface;
+class gp_Cone;
+class gp_Vec;
+class gp_XYZ;
 class IntPatch_Point;
 class IntSurf_PntOn2S;
 class math_Vector;
@@ -52,11 +55,17 @@ public:
   //! If theIsReversed is TRUE then theQSurf corresponds to the 
   //! second (otherwise, the first) surface while forming
   //! intersection point IntSurf_PntOn2S.
+  //! All math_Vector-objects must be filled as follows:
+  //!   [1] - U-parameter of thePSurf;
+  //!   [2] - V-parameter of thePSurf;
+  //!   [3] - U- (if V-isoline is considered) or V-parameter
+  //!         (if U-isoline is considered) of theQSurf.
   Standard_EXPORT static Standard_Boolean
                       AddPointOnUorVIso(const Handle(Adaptor3d_HSurface)& theQSurf,
                                         const Handle(Adaptor3d_HSurface)& thePSurf,
                                         const IntSurf_PntOn2S& theRefPt,
                                         const Standard_Boolean theIsU,
+                                        const Standard_Real theIsoParameter,
                                         const math_Vector& theToler,
                                         const math_Vector& theInitPoint,
                                         const math_Vector& theInfBound,
@@ -80,7 +89,6 @@ public:
                       AddSingularPole(const Handle(Adaptor3d_HSurface)& theQSurf,
                                       const Handle(Adaptor3d_HSurface)& thePSurf,
                                       const IntSurf_PntOn2S& thePtIso,
-                                      const Standard_Real theTol3d,
                                       IntPatch_Point& theVertex,
                                       IntSurf_PntOn2S& theAddedPoint,                                      
                                       const Standard_Boolean theIsReversed =
@@ -120,6 +128,39 @@ public:
                                            const Standard_Real theArrPeriods[4],
                                            IntSurf_PntOn2S &theNewPoint,
                                            IntPatch_Point* const theVertex = 0);
+
+protected:
+  //! Computes "special point" in the sphere
+  //! The parameter will be found in the range [0, 2*PI].
+  //! Therefore it must be adjusted to valid range by
+  //! the high-level algorithm
+  static Standard_EXPORT Standard_Boolean ProcessSphere(const IntSurf_PntOn2S& thePtIso,
+                                                        const gp_Vec& theDUofPSurf,
+                                                        const gp_Vec& theDVofPSurf,
+                                                        const Standard_Boolean theIsReversed,
+                                                        const Standard_Real theVquad,
+                                                        Standard_Real& theUquad,
+                                                        Standard_Boolean& theIsIsoChoosen);
+
+  //! Computes "special point" in the cone.
+  //! The parameter will be found in the range [0, 2*PI].
+  //! Therefore it must be adjusted to valid range by
+  //! the high-level algorithm.
+  static Standard_EXPORT Standard_Boolean ProcessCone(const IntSurf_PntOn2S& thePtIso,
+                                                      const gp_Vec& theDUofPSurf,
+                                                      const gp_Vec& theDVofPSurf,
+                                                      const gp_Cone& theCone,
+                                                      const Standard_Boolean theIsReversed,
+                                                      Standard_Real& theUquad,
+                                                      Standard_Boolean& theIsIsoChoosen);
+
+  //! Computes vector tangent to the intersection line in cone apex.
+  //! There exist not more than 2 tangent. They will be stores in theResult vector.
+  //! Returns the number of found tangents.
+  //! thePlnNormal is the normalized vector of the normal to the plane intersected the cone.
+  static Standard_EXPORT Standard_Integer GetTangentToIntLineForCone(const Standard_Real theConeSemiAngle,
+                                                                     const gp_XYZ& thePlnNormal,
+                                                                     gp_XYZ theResult[2]);
 };
 
 #endif // _IntPatch_AddSpecialPoints_HeaderFile

src/IntPatch/IntPatch_WLineTool.cxx

@@ -41,7 +41,8 @@ enum IntPatchWT_WLsConnectionType
 {
   IntPatchWT_NotConnected,
   IntPatchWT_Singular,
-  IntPatchWT_EachOther
+  IntPatchWT_Common,
+  IntPatchWT_ReqExtend
 };
 
 //=======================================================================
@@ -849,17 +850,23 @@ static IntPatchWT_WLsConnectionType
                                            const Standard_Real* const theArrPeriods)
 {
   const Standard_Real aSqToler = theToler3D*theToler3D;
-
+  IntPatchWT_WLsConnectionType aRetVal = IntPatchWT_NotConnected;
   if(theVec3.SquareMagnitude() <= aSqToler)
   {
-    return IntPatchWT_NotConnected;
+    if ((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle))
+    {
+      return aRetVal;
+    }
+    else
+    {
+      aRetVal = IntPatchWT_Common;
+    }
   }
-
+  else if((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle) ||
-  if((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle) ||
+          (theVec1.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle) ||
-     (theVec1.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle) ||
+          (theVec2.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle))
-     (theVec2.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle))
   {
-    return IntPatchWT_NotConnected;
+    return aRetVal;
   }
 
   const gp_Pnt aPmid(0.5*(thePtWL1.Value().XYZ()+thePtWL2.Value().XYZ()));
@@ -973,7 +980,12 @@ static IntPatchWT_WLsConnectionType
     return IntPatchWT_Singular;
   }
 
-  return IntPatchWT_EachOther;
+  if (aRetVal == IntPatchWT_Common)
+  {
+    return IntPatchWT_Common;
+  }
+
+  return IntPatchWT_ReqExtend;
 }
 
 //=======================================================================
@@ -984,12 +996,47 @@ static IntPatchWT_WLsConnectionType
 Standard_Boolean CheckArgumentsToJoin(const Handle(Adaptor3d_HSurface)& theS1,
                                       const Handle(Adaptor3d_HSurface)& theS2,
                                       const IntSurf_PntOn2S& thePnt,
+                                      const gp_Pnt& theP1,
+                                      const gp_Pnt& theP2,
+                                      const gp_Pnt& theP3,
                                       const Standard_Real theMinRad)
 {
   const Standard_Real aRad =
-        IntPatch_PointLine::CurvatureRadiusOfIntersLine(theS1, theS2, thePnt);
+    IntPatch_PointLine::CurvatureRadiusOfIntersLine(theS1, theS2, thePnt);
 
-  return (aRad > theMinRad);
+  if (aRad > theMinRad)
+  {
+    return Standard_True;
+  }
+  else if (aRad > 0.0)
+  {
+    return Standard_False;
+  }
+
+  // Curvature radius cannot be computed.
+  // Check smoothness of polygon.
+
+  //                  theP2
+  //                    *
+  //                    |
+  //                    |
+  //       *            o         *
+  //      theP1         O       theP3
+
+  //Joining is enabled if two conditions are satisfied together:
+  //  1. Angle (theP1, theP2, theP3) is quite big;
+  //  2. Modulus of perpendicular (O->theP2) to the segment (theP1->theP3)
+  //  is less than 0.01*<modulus of this segment>.
+
+  const gp_Vec aV12f(theP1, theP2), aV12l(theP2, theP3);
+
+  if (aV12f.Angle(aV12l) > IntPatch_WLineTool::myMaxConcatAngle)
+    return Standard_False;
+
+  const gp_Vec aV13(theP1, theP3);
+  const Standard_Real aSq13 = aV13.SquareMagnitude();
+
+  return (aV12f.CrossSquareMagnitude(aV13) < 1.0e-4*aSq13*aSq13);
 }
 
 //=======================================================================
@@ -1482,18 +1529,14 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       const Standard_Real aSqMinFDist = Min(aSqDistF, aSqDistL);
       if (aSqMinFDist < Precision::SquareConfusion())
       {
-        if (CheckArgumentsToJoin(theS1, theS2, aPntFWL1, aMinRad))
+        const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+        const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
+        const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) :
+                                             aWLine2->Point(aNbPntsWL2 - 1);
+        if (!IsSeamOrBound(aPt1, aPt2, aPntFWL1,
+                            anArrPeriods, anArrFBonds, anArrLBonds))
         {
-          const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+          isFirstConnected = Standard_True;
-          const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
-          const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) :
-                                                      aWLine2->Point(aNbPntsWL2 - 1);
-
-          if (!IsSeamOrBound(aPt1, aPt2, aPntFWL1,
-                              anArrPeriods, anArrFBonds, anArrLBonds))
-          {
-            isFirstConnected = Standard_True;
-          }
         }
       }
 
@@ -1503,18 +1546,14 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       const Standard_Real aSqMinLDist = Min(aSqDistF, aSqDistL);
       if (aSqMinLDist < Precision::SquareConfusion())
       {
-        if (CheckArgumentsToJoin(theS1, theS2, aPntLWL1, aMinRad))
+        const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+        const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1 - 1);
+        const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) :
+                                             aWLine2->Point(aNbPntsWL2 - 1);
+        if (!IsSeamOrBound(aPt1, aPt2, aPntLWL1,
+                           anArrPeriods, anArrFBonds, anArrLBonds))
         {
-          const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+          isLastConnected = Standard_True;
-          const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1 - 1);
-          const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) :
-                                               aWLine2->Point(aNbPntsWL2 - 1);
-
-          if (!IsSeamOrBound(aPt1, aPt2, aPntLWL1,
-                              anArrPeriods, anArrFBonds, anArrLBonds))
-          {
-            isLastConnected = Standard_True;
-          }
         }
       }
 
@@ -1549,15 +1588,29 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
 
     const Standard_Integer anIndexWL2 = isFirstConnected ? aListFC.First() : aListLC.First();
     Handle(IntPatch_WLine) aWLine2(Handle(IntPatch_WLine)::DownCast(theSlin.Value(anIndexWL2)));
-    
     const Standard_Integer aNbPntsWL2 = aWLine2->NbPnts();
     const IntSurf_PntOn2S& aPntFWL2 = aWLine2->Point(1);
-
+    const IntSurf_PntOn2S& aPntLWL2 = aWLine2->Point(aNbPntsWL2);
-    aWLine1->ClearVertexes();
     
     if (isFirstConnected)
     {
-      if (aPntFWL1.IsSame(aPntFWL2, Precision::Confusion()))
+      const Standard_Real aSqDistF = aPntFWL1.Value().SquareDistance(aPntFWL2.Value());
+      const Standard_Real aSqDistL = aPntFWL1.Value().SquareDistance(aPntLWL2.Value());
+      const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+
+      const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
+      const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) : 
+                                           aWLine2->Point(aNbPntsWL2 - 1);
+
+      if (!CheckArgumentsToJoin(theS1, theS2, aPntFWL1, aPt1.Value(),
+                                aPntFWL1.Value(), aPt2.Value(), aMinRad))
+      {
+        continue;
+      }
+
+      aWLine1->ClearVertexes();
+
+      if (isFM)
       {
         //First-First-connection
         for (Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
@@ -1578,10 +1631,26 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
     }
     else //if (isLastConnected)
     {
-      if (aPntLWL1.IsSame(aPntFWL2, Precision::Confusion()))
+      const Standard_Real aSqDistF = aPntLWL1.Value().SquareDistance(aPntFWL2.Value());
+      const Standard_Real aSqDistL = aPntLWL1.Value().SquareDistance(aPntLWL2.Value());
+
+      const Standard_Boolean isFM = (aSqDistF < aSqDistL);
+      const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1 - 1);
+      const IntSurf_PntOn2S& aPt2 = isFM ? aWLine2->Point(2) :
+                                           aWLine2->Point(aNbPntsWL2 - 1);
+
+      if (!CheckArgumentsToJoin(theS1, theS2, aPntLWL1, aPt1.Value(),
+                                aPntFWL1.Value(), aPt2.Value(), aMinRad))
+      {
+        continue;
+      }
+      
+      aWLine1->ClearVertexes();
+      
+      if (isFM)
       {
         //Last-First connection
-        for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
+        for (Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
         {
           const IntSurf_PntOn2S& aPt = aWLine2->Point(aNPt);
           aWLine1->Curve()->Add(aPt);
@@ -1609,13 +1678,15 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
 //purpose  : Performs extending theWLine1 and theWLine2 through their
 //            respecting end point.
 //=======================================================================
-void IntPatch_WLineTool::ExtendTwoWLines(IntPatch_SequenceOfLine& theSlin,
+void IntPatch_WLineTool::
-                                         const Handle(Adaptor3d_HSurface)& theS1,
+        ExtendTwoWLines(IntPatch_SequenceOfLine& theSlin,
-                                         const Handle(Adaptor3d_HSurface)& theS2,
+                        const Handle(Adaptor3d_HSurface)& theS1,
-                                         const Standard_Real theToler3D,
+                        const Handle(Adaptor3d_HSurface)& theS2,
-                                         const Standard_Real* const theArrPeriods,
+                        const Standard_Real theToler3D,
-                                         const Bnd_Box2d& theBoxS1,
+                        const Standard_Real* const theArrPeriods,
-                                         const Bnd_Box2d& theBoxS2)
+                        const Bnd_Box2d& theBoxS1,
+                        const Bnd_Box2d& theBoxS2,
+                        const NCollection_List<gp_Pnt>& theListOfCriticalPoints)
 {
   if(theSlin.Length() < 2)
     return;
@@ -1677,6 +1748,46 @@ void IntPatch_WLineTool::ExtendTwoWLines(IntPatch_SequenceOfLine& theSlin,
       {
         aCheckResult |= IntPatchWT_DisLastFirst | IntPatchWT_DisLastLast;
       }
+
+      if (!theListOfCriticalPoints.IsEmpty())
+      {
+        for (NCollection_List<gp_Pnt>::Iterator anItr(theListOfCriticalPoints);
+             anItr.More(); anItr.Next())
+        {
+          const gp_Pnt &aPt = anItr.Value();
+          if (!(aCheckResult & (IntPatchWT_DisFirstFirst | IntPatchWT_DisFirstLast)))
+          {
+            if (aPt.SquareDistance(aPntFWL1.Value()) < Precision::Confusion())
+            {
+              aCheckResult |= IntPatchWT_DisFirstFirst | IntPatchWT_DisFirstLast;
+            }
+          }
+
+          if (!(aCheckResult & (IntPatchWT_DisLastFirst | IntPatchWT_DisLastLast)))
+          {
+            if (aPt.SquareDistance(aPntLWL1.Value()) < Precision::Confusion())
+            {
+              aCheckResult |= IntPatchWT_DisLastFirst | IntPatchWT_DisLastLast;
+            }
+          }
+
+          if (!(aCheckResult & (IntPatchWT_DisFirstFirst | IntPatchWT_DisLastFirst)))
+          {
+            if (aPt.SquareDistance(aPntFWL2.Value()) < Precision::Confusion())
+            {
+              aCheckResult |= IntPatchWT_DisFirstFirst | IntPatchWT_DisLastFirst;
+            }
+          }
+
+          if (!(aCheckResult & (IntPatchWT_DisFirstLast | IntPatchWT_DisLastLast)))
+          {
+            if (aPt.SquareDistance(aPntLWL2.Value()) < Precision::Confusion())
+            {
+              aCheckResult |= IntPatchWT_DisFirstLast | IntPatchWT_DisLastLast;
+            }
+          }
+        }
+      }
     }
 
     if(aCheckResult == (IntPatchWT_DisFirstFirst | IntPatchWT_DisFirstLast |

src/IntPatch/IntPatch_WLineTool.hxx

@@ -18,6 +18,7 @@
 
 #include <IntPatch_SequenceOfLine.hxx>
 #include <IntPatch_WLine.hxx>
+#include <NCollection_List.hxx>
 
 class Adaptor3d_TopolTool;
 
@@ -49,29 +50,30 @@ public:
                                               const Handle(Adaptor3d_TopolTool)  &theDom1,
                                               const Handle(Adaptor3d_TopolTool)  &theDom2);
 
-//! Joins all WLines from theSlin to one if it is possible and records 
+  //! Joins all WLines from theSlin to one if it is possible and records 
-//! the result into theSlin again. Lines will be kept to be splitted if:
+  //! the result into theSlin again. Lines will be kept to be split if:
-//! a) they are separated (has no common points);
+  //! a) they are separated (has no common points);
-//! b) resulted line (after joining) go through seam-edges or surface boundaries.
+  //! b) resulted line (after joining) go through seam-edges or surface boundaries.
-//!
+  //!
-//! In addition, if points in theSPnt lies at least in one of the line in theSlin,
+  //! In addition, if points in theSPnt lies at least in one of the line in theSlin,
-//! this point will be deleted.
+  //! this point will be deleted.
   Standard_EXPORT static void JoinWLines(IntPatch_SequenceOfLine& theSlin,
                                          IntPatch_SequenceOfPoint& theSPnt,
                                          Handle(Adaptor3d_HSurface) theS1,
                                          Handle(Adaptor3d_HSurface) theS2,
                                          const Standard_Real theTol3D);
 
-//! Extends every line from theSlin (if it is possible) to be started/finished
+  //! Extends every line from theSlin (if it is possible) to be started/finished
-//! in strictly determined point (in the place of joint of two lines).
+  //! in strictly determined point (in the place of joint of two lines).
-//! As result, some gaps between two lines will vanish.
+  //! As result, some gaps between two lines will vanish.
-//! The Walking lines are supposed (algorithm will do nothing for not-Walking line)
+  //! The Walking lines are supposed (algorithm will do nothing for not-Walking line)
-//! to be computed as a result of intersection. Both theS1 and theS2 
+  //! to be computed as a result of intersection. Both theS1 and theS2 
-//! must be quadrics. Other cases are not supported.
+  //! must be quadrics. Other cases are not supported.
-//! theArrPeriods must be filled as follows (every value must not be negative;
+  //! theArrPeriods must be filled as follows (every value must not be negative;
-//! if the surface is not periodic the period must be equal to 0.0 strictly):
+  //! if the surface is not periodic the period must be equal to 0.0 strictly):
-//! {<U-period of 1st surface>, <V-period of 1st surface>,
+  //! {<U-period of 1st surface>, <V-period of 1st surface>,
-//!               <U-period of 2nd surface>, <V-period of 2nd surface>}.
+  //!               <U-period of 2nd surface>, <V-period of 2nd surface>}.
+  //! theListOfCriticalPoints must contain 3D-points where joining is disabled.
   Standard_EXPORT static void
             ExtendTwoWLines(IntPatch_SequenceOfLine& theSlin,
                             const Handle(Adaptor3d_HSurface)& theS1,
@@ -79,7 +81,8 @@ public:
                             const Standard_Real theToler3D,
                             const Standard_Real* const theArrPeriods,
                             const Bnd_Box2d& theBoxS1,
-                            const Bnd_Box2d& theBoxS2);
+                            const Bnd_Box2d& theBoxS2,
+                            const NCollection_List<gp_Pnt>& theListOfCriticalPoints);
 
   //! Max angle to concatenate two WLines to avoid result with C0-continuity
   static const Standard_Real myMaxConcatAngle;

src/IntStart/IntStart_SearchOnBoundaries.gxx

@@ -225,14 +225,13 @@ void BoundedArc (const TheArc& A,
   // Creer l echantillonage (math_FunctionSample ou classe heritant)
   // Appel a math_FunctionAllRoots
 
-  Standard_Real EpsX = TheArcTool::Resolution(A,Precision::Confusion());
   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
   //@@@ La Tolerance est asociee a l arc  ( Incoherence avec le cheminement )
   //@@@   ( EpsX ~ 1e-5   et ResolutionU et V ~ 1e-9 )
   //@@@   le vertex trouve ici n'est pas retrouve comme point d arret d une 
   //@@@   ligne de cheminement
   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
-  EpsX = 0.0000000001;
+  Standard_Real EpsX = 1.e-10;
   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

src/QABugs/QABugs_19.cxx

@@ -1533,25 +1533,6 @@ static Standard_Integer OCC24271 (Draw_Interpretor& di,
 #include <GeomInt_IntSS.hxx>
 #include <Geom_ConicalSurface.hxx>
 #include <Standard_ErrorHandler.hxx>
-//=======================================================================
-//function : OCC23972
-//purpose  : 
-//=======================================================================
-static void DoGeomIntSSTest (const Handle(Geom_Surface)& theSurf1,
-			     const Handle(Geom_Surface)& theSurf2,
-			     const Standard_Integer theNbSol,
-			     Draw_Interpretor& di)
-{
-  try {
-    OCC_CATCH_SIGNALS
-	 GeomInt_IntSS anInter;
-	 anInter.Perform (theSurf1, theSurf2, Precision::Confusion(), Standard_True);
-	 QVERIFY (anInter.IsDone());
-	 QCOMPARE (anInter.NbLines(), theNbSol);
-  } catch (...) {
-    QVERIFY (Standard_False);
-  }
-}
 
 namespace {
   static Handle(Geom_ConicalSurface) CreateCone (const gp_Pnt& theLoc,
@@ -1568,28 +1549,30 @@ namespace {
   }
 }
 
-static Standard_Integer OCC23972 (Draw_Interpretor& di,Standard_Integer n, const char**)
+static Standard_Integer OCC23972(Draw_Interpretor& /*theDI*/,
+                                 Standard_Integer theNArg, const char** theArgs)
 {
-  if (n != 1) return 1;
+  if (theNArg != 3) return 1;
 
-  //process specific cones, cannot read them from files because due to rounding the original error
+  //process specific cones, cannot read them from files because 
-  //in math_FunctionRoots gets hidden
+  //due to rounding the original error in math_FunctionRoots gets hidden
-  Handle(Geom_Surface) aS1 = CreateCone (
+  const Handle(Geom_Surface) aS1 = CreateCone(
-					 gp_Pnt (123.694345356663, 789.9, 68.15),
+                              gp_Pnt(123.694345356663, 789.9, 68.15),
-					 gp_Dir (-1, 3.48029791472957e-016, -8.41302743359754e-017),
+                              gp_Dir(-1, 3.48029791472957e-016, -8.41302743359754e-017),
-					 gp_Dir (-3.48029791472957e-016, -1, -3.17572289932207e-016),
+                              gp_Dir(-3.48029791472957e-016, -1, -3.17572289932207e-016),
-					 3.28206830417112,
+                              3.28206830417112,
-					 0.780868809443031,
+                              0.780868809443031,
-					 0.624695047554424);
+                              0.624695047554424);
-  Handle(Geom_Surface) aS2 = CreateCone (
+  const Handle(Geom_Surface) aS2 = CreateCone(
-					 gp_Pnt (123.694345356663, 784.9, 68.15),
+                              gp_Pnt(123.694345356663, 784.9, 68.15),
-					 gp_Dir (-1, -2.5209507537117e-016, -1.49772808948866e-016),
+                              gp_Dir(-1, -2.5209507537117e-016, -1.49772808948866e-016),
-					 gp_Dir (1.49772808948866e-016, 3.17572289932207e-016, -1),
+                              gp_Dir(1.49772808948866e-016, 3.17572289932207e-016, -1),
-					 3.28206830417112,
+                              3.28206830417112,
-					 0.780868809443031,
+                              0.780868809443031,
-					 0.624695047554424);
+                              0.624695047554424);
   
-  DoGeomIntSSTest (aS1, aS2, 2, di);
+  DrawTrSurf::Set(theArgs[1], aS1);
+  DrawTrSurf::Set(theArgs[2], aS2);
 
   return 0;
 }

src/QABugs/QABugs_20.cxx

@@ -2891,6 +2891,56 @@ static Standard_Integer OCC29531(Draw_Interpretor&, Standard_Integer, const char
   return 0;
 }
 
+//=======================================================================
+//function : OCC29807
+//purpose  : 
+//=======================================================================
+#include <GeomAdaptor_HSurface.hxx>
+#include <IntPatch_PointLine.hxx>
+#include <IntSurf_PntOn2S.hxx>
+static Standard_Integer OCC29807(Draw_Interpretor& theDI, Standard_Integer theNArg, const char** theArgV)
+{
+  if (theNArg != 7)
+  {
+    theDI << "Use: " << theArgV[0] << "surface1 surface2 u1 v1 u2 v2\n";
+    return 1;
+  }
+
+  const Handle(Geom_Surface) aS1 = DrawTrSurf::GetSurface(theArgV[1]);
+  const Handle(Geom_Surface) aS2 = DrawTrSurf::GetSurface(theArgV[2]);
+
+  if (aS1.IsNull() || aS2.IsNull())
+  {
+    theDI << "Error. Null surface is not supported.\n";
+    return 1;
+  }
+
+  const Standard_Real aU1 = Draw::Atof(theArgV[3]);
+  const Standard_Real aV1 = Draw::Atof(theArgV[4]);
+  const Standard_Real aU2 = Draw::Atof(theArgV[5]);
+  const Standard_Real aV2 = Draw::Atof(theArgV[6]);
+
+  const Handle(GeomAdaptor_HSurface) anAS1 = new GeomAdaptor_HSurface(aS1);
+  const Handle(GeomAdaptor_HSurface) anAS2 = new GeomAdaptor_HSurface(aS2);
+
+  const gp_Pnt aP1 = anAS1->Value(aU1, aV1);
+  const gp_Pnt aP2 = anAS2->Value(aU2, aV2);
+
+  if (aP1.SquareDistance(aP2) > Precision::SquareConfusion())
+  {
+    theDI << "Error. True intersection point must be specified. "
+             "Please check parameters: u1 v1 u2 v2.\n";
+    return 1;
+  }
+
+  IntSurf_PntOn2S aPOn2S;
+  aPOn2S.SetValue(0.5*(aP1.XYZ() + aP2.XYZ()), aU1, aV1, aU2, aV2);
+
+  const Standard_Real aCurvatureRadius = IntPatch_PointLine::CurvatureRadiusOfIntersLine(anAS1, anAS2, aPOn2S);
+  theDI << "Radius of curvature is " << aCurvatureRadius << "\n";
+  return 0;
+}
+
 //=======================================================================
 //function : OCC29925
 //purpose  : check safety of functions like IsSpace(), LowerCase(), etc. for all chars
@@ -2977,5 +3027,7 @@ void QABugs::Commands_20(Draw_Interpretor& theCommands) {
 
   theCommands.Add ("OCC29064", "OCC29064: test memory usage by copying empty maps", __FILE__, OCC29064, group);
   theCommands.Add ("OCC29925", "OCC29925: check safety of character classification functions", __FILE__, OCC29925, group);
+  theCommands.Add("OCC29807", "OCC29807 surface1 surface2 u1 v1 u2 v2", __FILE__, OCC29807, group);
+
   return;
 }

src/gp/gp_Cone.hxx

@@ -35,7 +35,7 @@ class gp_Vec;
 
 
 //! Defines an infinite conical surface.
-//! A cone is defined by its half-angle at the apex and
+//! A cone is defined by its half-angle (can be negative) at the apex and
 //! positioned in space with a coordinate system (a gp_Ax3
 //! object) and a "reference radius" where:
 //! -   the "main Axis" of the coordinate system is the axis of   revolution of the cone,
@@ -76,12 +76,13 @@ public:
 
   //! Creates an infinite conical surface. A3 locates the cone
   //! in the space and defines the reference plane of the surface.
-  //! Ang is the conical surface semi-angle between 0 and PI/2 radians.
+  //! Ang is the conical surface semi-angle. Its absolute value is in range
+  //! ]0, PI/2[.
   //! Radius is the radius of the circle in the reference plane of
   //! the cone.
   //! Raises ConstructionError
-  //! . if Radius is lower than 0.0
+  //! * if Radius is lower than 0.0
-  //! . Ang < Resolution from gp  or Ang >= (PI/2) - Resolution.
+  //! * Abs(Ang) < Resolution from gp  or Abs(Ang) >= (PI/2) - Resolution.
     gp_Cone(const gp_Ax3& A3, const Standard_Real Ang, const Standard_Real Radius);
   
   //! Changes the symmetry axis of the cone.  Raises ConstructionError
@@ -105,8 +106,9 @@ public:
   
 
   //! Changes the semi-angle of the cone.
-  //! Ang is the conical surface semi-angle  ]0,PI/2[.
+  //! Semi-angle can be negative. Its absolute value
-  //! Raises ConstructionError if Ang < Resolution from gp or Ang >= PI/2 - Resolution
+  //! Abs(Ang) is in range ]0,PI/2[.
+  //! Raises ConstructionError if Abs(Ang) < Resolution from gp or Abs(Ang) >= PI/2 - Resolution
     void SetSemiAngle (const Standard_Real Ang);
   
 
@@ -146,6 +148,7 @@ public:
     Standard_Real RefRadius() const;
   
   //! Returns the half-angle at the apex of this cone.
+  //! Attention! Semi-angle can be negative.
     Standard_Real SemiAngle() const;
   
   //! Returns the XAxis of the reference plane.

tests/boolean/volumemaker/C5

@@ -1,10 +1,6 @@
 # test script on make volume operation
 # cone plane
 
-puts "TODO CR28503 ALL: Error : The area of result shape is"
-puts "TODO CR28503 ALL: Error :  is WRONG because number of SOLID entities in shape"
-puts "TODO CR28503 ALL: Faulty shapes in variables faulty_"
-
 # planar face 
 plane pln_f1 27.577164466275352 -1038.2137499999999 27.577164466275359 0.70710678118654746 4.4408920985006262e-016 0.70710678118654768
 erase pln_f1
@@ -38,6 +34,6 @@ mkface f6 con_f6 0 6.2831853071795862 0 1000000
 # make volume operation 
 mkvolume result f1 f2 f3 f4 f5 f6
 
-checkprops result -s 5.19571e+006
+checkprops result -s 5.5768e+006
 
-checknbshapes result -solid 16
+checknbshapes result -solid 17

tests/boolean/volumemaker/C6

@@ -1,8 +1,6 @@
 # test script on make volume operation
 # cone plane
 
-puts "TODO CR28503 ALL: Error :  is WRONG because number of SOLID entities in shape"
-
 # planar face 
 plane pln_f1 -306.53078964627537 -1038.2137499999999 -251.37646071372467 -0.70710678118654746 4.4408920985006262e-016 0.70710678118654768
 erase pln_f1
@@ -43,4 +41,4 @@ mkvolume result f1 f2 f3 f4 f5 f6 f7
 
 checkprops result -s 6.45353e+006
 
-checknbshapes result -solid 29
+checknbshapes result -solid 30

tests/boolean/volumemaker/E7

@@ -29,5 +29,6 @@ mkface f5 con_f5 0 6.2831853071795862 0 1000000
 # make volume operation 
 mkvolume result f1 f2 f3 f4 f5
 
-checkprops result -s 5.28759e+006
+checkprops result -s 5.7053e+006
 
+checknbshapes result -solid 8

tests/bugs/fclasses/bug23972

@@ -1,11 +0,0 @@
-puts "=========="
-puts "OCC23972"
-puts "=========="
-puts ""
-###########################################################################
-## Exception thrown when intersecting two cones
-###########################################################################
-
-pload QAcommands
-
-OCC23972

tests/bugs/modalg_2/bug20964_1

@@ -1,29 +1,38 @@
 puts "============"
-puts "OCC20964"
+puts "OCC20964: Wrong result of cut operation for given shapes"
 puts "============"
 puts ""
-#######################################################################
-# Wrong result of cut operation for given shapes
-#######################################################################
 
-set BugNumber OCC20964
-
-puts "Load first shape ..."
 restore [locate_data_file OCC20964_revsolid.brep] b1 
-puts "Load second shape ..."
 restore [locate_data_file OCC20964_sphere.brep] b2 
 
-puts "Prepare boolean operation ..."
+bclearobjects
-bop b1 b2
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
 
-puts "Start boolean operation ..."
+# SECTION
-bopsection result
+bbop result 4
-puts "Finish boolean operation ..."
 
-checkprops result -l 323.636 
+checkprops result -l 323.635
 checkshape result 
-checksection result
+checksection result -r 0
 
-checknbshapes result -vertex 6 -edge 6 -wire 0 -face 0 -shell 0 -solid 0 -compsolid 0 -compound 1 -shape 13
+set NbShapesRef "
+Number of shapes in .*
+ VERTEX : 6
+ EDGE : 6
+ WIRE : 0
+ FACE : 0
+ SHELL : 0
+ SOLID : 0
+ COMPSOLID : 0
+ COMPOUND : 1
+ SHAPE : 13
+"
+
+checknbshapes result -ref $NbShapesRef
+checkmaxtol result -ref 0.013928665225777443
 
 checkview -display result -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_2/bug20964_2

@@ -1,28 +1,48 @@
 puts "============"
-puts "OCC20964"
+puts "OCC20964: Wrong result of cut operation for given shapes"
 puts "============"
 puts ""
-#######################################################################
-# Wrong result of cut operation for given shapes
-#######################################################################
 
-set BugNumber OCC20964
-
-puts "Load first shape ..."
 restore [locate_data_file OCC20964_revsolid.brep] b1 
-puts "Load second shape ..."
 restore [locate_data_file OCC20964_sphere.brep] b2 
 
-puts "Prepare boolean operation ..."
+bclearobjects
-bop b1 b2
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
 
-puts "Start boolean operation ..."
+# COMMON
-bopcommon result
+bbop result 0
-puts "Finish boolean operation ..."
 
-checkprops result -s 5164.66 
+checkprops result -s 5158.93
 checkshape result
 
-checknbshapes result -vertex 8 -edge 12 -wire 6 -face 5 -shell 1 -solid 1 -compsolid 0 -compound 1 -shape 34
+set NbShapesRef "
+Number of shapes in .*
+ VERTEX : 8
+ EDGE : 12
+ WIRE : 6
+ FACE : 5
+ SHELL : 1
+ SOLID : 1
+ COMPSOLID : 0
+ COMPOUND : 1
+ SHAPE : 34
+"
 
-checkview -display result -2d -path ${imagedir}/${test_image}.png
+checknbshapes result -ref $NbShapesRef
+
+checkmaxtol result -ref 0.013928665225777443
+
+checkview -display result -2d -path ${imagedir}/${test_image}_axo.png
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 15.9695 -proj 0.0161244 -0.535424 0.84443 -up -0.773936 0.528029 0.349583 -at 117.532 248.227 -4.41145e-007 -eye 118.783 206.661 65.5549
+
+if { [string compare "" [tricheck result] ] } {
+  puts "Error in triangulation"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}_shade.png

tests/bugs/modalg_2/bug20964_3

@@ -1,28 +1,48 @@
 puts "============"
-puts "OCC20964"
+puts "OCC20964: Wrong result of cut operation for given shapes"
 puts "============"
 puts ""
-#######################################################################
-# Wrong result of cut operation for given shapes
-#######################################################################
 
-set BugNumber OCC20964
-
-puts "Load first shape ..."
 restore [locate_data_file OCC20964_revsolid.brep] b1 
-puts "Load second shape ..."
 restore [locate_data_file OCC20964_sphere.brep] b2 
 
-puts "Prepare boolean operation ..."
+bclearobjects
-bop b1 b2
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
 
-puts "Start boolean operation ..."
+# FUSE
-bopfuse result
+bbop result 1
-puts "Finish boolean operation ..."
 
-checkprops result -s 828829 
+checkprops result -s 828808
 checkshape result
 
-checknbshapes result -vertex 20 -edge 32 -wire 18 -face 16 -shell 3 -solid 1 -compsolid 0 -compound 1 -shape 91
+set NbShapesRef "
+Number of shapes in .*
+ VERTEX : 20
+ EDGE : 32
+ WIRE : 18
+ FACE : 16
+ SHELL : 3
+ SOLID : 1
+ COMPSOLID : 0
+ COMPOUND : 1
+ SHAPE : 91
+"
 
-checkview -display result -2d -path ${imagedir}/${test_image}.png
+checknbshapes result -ref $NbShapesRef
+
+checkmaxtol result -ref 0.013928665225777443
+
+checkview -display result -2d -path ${imagedir}/${test_image}_axo.png
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 5.62548 -proj 0.751453 0.00874872 0.659728 -up -0.600015 0.42493 0.677803 -at 99.378 322.007 8.81412 -eye 467.31 326.291 331.835
+
+if { [string compare "" [tricheck result] ] } {
+  puts "Error in triangulation"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}_shade.png

tests/bugs/modalg_2/bug20964_4

@@ -1,28 +1,47 @@
 puts "============"
-puts "OCC20964"
+puts "OCC20964: Wrong result of cut operation for given shapes"
 puts "============"
 puts ""
-#######################################################################
-# Wrong result of cut operation for given shapes
-#######################################################################
 
-set BugNumber OCC20964
-
-puts "Load first shape ..."
 restore [locate_data_file OCC20964_revsolid.brep] b1 
-puts "Load second shape ..."
 restore [locate_data_file OCC20964_sphere.brep] b2 
 
-puts "Prepare boolean operation ..."
+bclearobjects
-bop b1 b2
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
 
-puts "Start boolean operation ..."
+# CUT 1-2
-bopcut result
+bbop result 2
-puts "Finish boolean operation ..."
 
-checkprops result -s 821892 
+checkprops result -s 821864
 checkshape result
 
-checknbshapes result -vertex 22 -edge 35 -wire 17 -face 14 -shell 2 -solid 1 -compsolid 0 -compound 1 -shape 92
+set NbShapesRef "
+Number of shapes in .*
+ VERTEX : 22
+ EDGE : 35
+ WIRE : 17
+ FACE : 14
+ SHELL : 2
+ SOLID : 1
+ COMPSOLID : 0
+ COMPOUND : 1
+ SHAPE : 92
+"
 
-checkview -display result -2d -path ${imagedir}/${test_image}.png
+checknbshapes result -ref $NbShapesRef
+checkmaxtol result -ref 0.013928665225777443
+
+checkview -display result -2d -path ${imagedir}/${test_image}_axo.png
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 3.80971 -proj 0.769527 -0.23844 0.59243 -up -0.456402 0.443531 0.771348 -at 70.5761 168.919 28.0572 -eye 455.34 49.6989 324.272
+
+if { [string compare "" [tricheck result] ] } {
+  puts "Error in triangulation"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}_shade.png

tests/bugs/modalg_2/bug20964_5

@@ -1,28 +1,52 @@
 puts "============"
-puts "OCC20964"
+puts "OCC20964 Wrong result of cut operation for given shapes"
 puts "============"
 puts ""
-#######################################################################
-# Wrong result of cut operation for given shapes
-#######################################################################
 
-set BugNumber OCC20964
-
-puts "Load first shape ..."
 restore [locate_data_file OCC20964_revsolid.brep] b1 
-puts "Load second shape ..."
 restore [locate_data_file OCC20964_sphere.brep] b2 
 
-puts "Prepare boolean operation ..."
+bclearobjects
-bop b1 b2
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
 
-puts "Start boolean operation ..."
+# CUT 2-1
-boptuc result
+bbop result 3
-puts "Finish boolean operation ..."
 
-checkprops result -s 12101.2 
+checkprops result -s 12102.9
 checkshape result
 
-checknbshapes result -vertex 6 -edge 9 -wire 7 -face 7 -shell 2 -solid 2 -compsolid 0 -compound 1 -shape 34
+set NbShapesRef "
+Number of shapes in .*
+ VERTEX : 6
+ EDGE : 9
+ WIRE : 7
+ FACE : 7
+ SHELL : 2
+ SOLID : 2
+ COMPSOLID : 0
+ COMPOUND : 1
+ SHAPE : 34
+"
 
-checkview -display result -2d -path ${imagedir}/${test_image}.png
+checknbshapes result -ref $NbShapesRef
+
+checkmaxtol result -ref 0.013928665225777443
+
+checkview -display result -2d -path ${imagedir}/${test_image}_axo.png
+
+vdisplay result
+vsetdispmode 1
+vsetcolor result red
+vdisplay b1
+vsettransparency b1 0.5
+
+vviewparams -scale 11.5636 -proj 0.296876 -0.306332 0.904447 -up -0.63687 0.642216 0.426562 -at 120.352 250.434 3.97104e-006 -eye 147.307 222.621 82.1187
+
+if { [string compare "" [tricheck result] ] } {
+  puts "Error in triangulation"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}_shade.png

tests/bugs/modalg_6/bug27269

@@ -10,6 +10,10 @@ puts ""
 restore [locate_data_file bug27267_cmpd.brep] a
 explode a f
 
+smallview
+don a_7; fit
+disp a_1
+
 #############################
 set log [bopcurves a_1 a_7 -2d]
 #############################
@@ -33,4 +37,4 @@ for {set i 1} {$i <= ${NbCurv}} {incr i} {
   }
 }
 
-
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_6/bug27282_2

@@ -6,7 +6,7 @@ puts ""
 ## [Regression to 6.9.1] smesh/bugs_00/A6: Cut produces an empty shape
 ###############################
 
-set MaxTol 2.1243683206633536e-006
+set MaxTol 2.9376013151287501e-006
 set GoodNbCurv 1
 
 restore [locate_data_file bug27282_cmpd.brep] a

tests/bugs/modalg_6/bug28626_1

@@ -18,3 +18,8 @@ checknbshapes result -solid 1 -shell 1 -face 10 -wire 10 -edge 20 -vertex 12
 checkmaxtol result -ref 8.00001e-007
 
 checkprops result -v 1.35999e+006
+
+smallview
+don b result
+fit
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_6/bug28626_2

@@ -12,8 +12,6 @@ trimv tc1 c1 0 42.4264068711929
 trimv tc2 c2 0 42.4264068711929
 mkface f1 tc1
 mkface f2 tc2
-donly f1 f2
-fit
 
 bop f1 f2
 bopsection result
@@ -25,3 +23,8 @@ checknbshapes result -edge 2 -vertex 3
 checkmaxtol result -ref 5.21731e-007
 
 checkprops result -l 88.9692
+
+smallview
+don f1 f2 result
+fit
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_6/bug28626_3

@@ -25,3 +25,8 @@ checknbshapes result -edge 5 -vertex 6
 checkmaxtol result -ref 6.02982e-007
 
 checkprops result -l 94.3164
+
+smallview
+don f1 f2 result
+fit
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug23176

@@ -11,33 +11,51 @@ restore [locate_data_file bug23176_surface2_draw.draw] s2
 
 intersect result s1 s2
 
-# first result curve
+set che [whatis result]
-cvalue result_1 0 x y z
+set ind [string first "3d curve" $che]
-vertex v0 x y z
+if {${ind} >= 0} {
-cvalue result_1 1 x y z
+  #Only variable "result" exists
-vertex v1 x y z
+  renamevar result result_1
+}
 
-# second result curve
+bclearobjects
-cvalue result_2 0 x y z
+bcleartools
-vertex v2 x y z
-cvalue result_2 1 x y z
-vertex v3 x y z
 
-set tol_abs_dist 1.0e-7
+set ic 1
-set tol_rel_dist 0.001
+set AllowRepeat 1
+while { $AllowRepeat != 0 } {
+  set che [whatis result_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeat 0
+  } else {
    
-distmini d_f v0 v1
+    bounds result_$ic U1 U2
-regexp {([-0-9.+eE]+)} [dump d_f_val] full dist_1
-set expected_dist_1 408.093320004435
-checkreal "Length first curve" ${dist_1} ${expected_dist_1} ${tol_abs_dist} ${tol_rel_dist}
     
-distmini d_s v2 v3
+    if {[dval U2-U1] < 1.0e-20} {
-regexp {([-0-9.+eE]+)} [dump d_s_val] full dist_2
+      puts "Error: Wrong curve's range!"
-set expected_dist_2 408.093320004435
+    }
-checkreal "Length second curve" ${dist_2} ${expected_dist_2} ${tol_abs_dist} ${tol_rel_dist}
+    
+    xdistcs result_$ic s1 U1 U2 10 2.0e-7
+    xdistcs result_$ic s2 U1 U2 10 2.0e-7
+    
+    mkedge ee result_$ic
+    baddobjects ee
+    incr ic
+  }
+}
+
+# Check gaps in result
+bfillds
+bbuild rw
+
+checksection rw -r 0
+checkmaxtol rw -min_tol 2.0e-7
+
+checknbshapes rw -edge 4 -vertex 4
 
 smallview
-donly result*
+don result*
 fit
 display s1 s2
 checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug23972

@@ -0,0 +1,52 @@
+puts "=========="
+puts "OCC23972: Excep-tion thrown when intersecting two cones"
+puts "=========="
+puts ""
+
+pload QAcommands
+
+set GoodNbCurv 2
+
+OCC23972 s1 s2
+
+intersect res s1 s2
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  renamevar res res_1
+}
+
+set ic 1
+set AllowRepeat 1
+while { $AllowRepeat != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeat 0
+  } else {
+   
+    bounds res_$ic U1 U2
+    
+    # Hyperbola is expected as intersection result.
+    # So, all bounds are adjusted to allow correct computation.
+    if { [dval U1] < -20.0 } { dset U1 -20.0 }
+    if { [dval U2] > 20.0 }  { dset U2 20.0 }
+    
+    if {[dval U2-U1] < 1.0e-20} {
+      puts "Error: Wrong curve's range!"
+    }
+    
+    xdistcs res_$ic s1 U1 U2 10 3.0e-7
+    xdistcs res_$ic s2 U1 U2 10 3.0e-7
+    
+    incr ic
+  }
+}
+
+if {[expr {$ic - 1}] == $GoodNbCurv} {
+  puts "OK: Number of curves is good!"
+} else {
+  puts "Error: Number of curves is bad!"
+}

tests/bugs/modalg_7/bug25542

@@ -0,0 +1,111 @@
+puts "============"
+puts "OCC25542: Boolean operation fai-lure for Cylinder and Cone in critical location."
+puts "============"
+puts ""
+
+restore [locate_data_file bug25542_cylinder.brep] b1
+restore [locate_data_file bug25542_cone.brep] b2
+
+smallview
+
+bop b1 b2
+
+puts ""
+puts "FUSE"
+
+bopfuse rfu
+donly rfu
+fit
+xwd ${imagedir}/${casename}_rfu.png
+
+puts ""
+puts "COMMON"
+
+bopcommon rco
+donly rco
+fit
+xwd ${imagedir}/${casename}_rco.png
+
+puts ""
+puts "CUT"
+
+bopcut rcu
+donly rcu
+fit
+xwd ${imagedir}/${casename}_rcu.png
+
+puts ""
+puts "TUC"
+
+boptuc rtu
+donly rtu
+fit
+xwd ${imagedir}/${casename}_rtu.png
+
+
+set exp_nbshapes_rfu "
+Number of shapes in shape
+ VERTEX    : 6
+ EDGE      : 14
+ WIRE      : 8
+ FACE      : 8
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 39
+"
+
+set exp_nbshapes_rco "
+Number of shapes in shape
+ VERTEX    : 3
+ EDGE      : 5
+ WIRE      : 3
+ FACE      : 3
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 17
+"
+
+set exp_nbshapes_rcu "
+Number of shapes in shape
+ VERTEX    : 5
+ EDGE      : 9
+ WIRE      : 5
+ FACE      : 5
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 27
+"
+
+set exp_nbshapes_rtu "
+Number of shapes in shape
+ VERTEX    : 4
+ EDGE      : 8
+ WIRE      : 4
+ FACE      : 4
+ SHELL     : 1
+ SOLID     : 1
+ COMPSOLID : 0
+ COMPOUND  : 1
+ SHAPE     : 23
+"
+
+checknbshapes rfu -ref ${exp_nbshapes_rfu} -t -m "FUSE"
+checknbshapes rco -ref ${exp_nbshapes_rco} -t -m "COMMON"
+checknbshapes rcu -ref ${exp_nbshapes_rcu} -t -m "CUT"
+checknbshapes rtu -ref ${exp_nbshapes_rtu} -t -m "TUC"
+
+checkshape rfu
+checkshape rco
+checkshape rcu
+checkshape rtu
+
+checkprops rfu -s 59099.9
+checkprops rco -s 6951.3
+checkprops rcu -s 57145.3
+checkprops rtu -s 7759.96

tests/bugs/modalg_7/bug29807_b1

@@ -0,0 +1,22 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9465.07 -s 4012.74
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b2

@@ -0,0 +1,25 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+removeloc b2 b2
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9463.99 -s 4014.54
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b3a

@@ -0,0 +1,55 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+puts "TODO OCC29922 ALL: Error: Degenerated edge is not found"
+puts "TODO OCC29922 ALL: Error: Result contains not triangulated face"
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc b2 b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
+
+# CUT operation
+bbop result 2
+
+savehistory ResHist
+
+checkshape result
+checkprops result -v 9465.07 -s 4012.74
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+explode b2 f
+modified reshm ResHist b2_2
+
+checknbshapes reshm -face 2
+
+foreach a [explode reshm f] {
+  # Every modified face must contain degenerated edge
+  
+  if { ![regexp {degenerated} [dump $a] ] } {
+    puts "Error: Degenerated edge is not found"
+  }
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+if { [regexp {no triangulation} [tricheck result] ] } {
+  puts "Error: Result contains not triangulated face"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b3b

@@ -0,0 +1,57 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc b2 b2
+
+explode b2 v
+settolerance b2_2 1.0e-7
+
+checkshape b2
+
+bclearobjects
+bcleartools
+baddobjects b1
+baddtools b2
+bfillds
+
+# CUT operation
+bbop result 2
+
+savehistory ResHist
+
+checkshape result
+checkprops result -v 9465.07 -s 4012.74
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+explode b2 f
+modified reshm ResHist b2_2
+
+checknbshapes reshm -face 2
+
+foreach a [explode reshm f] {
+  # Every modified face must contain degenerated edge
+  
+  if { ![regexp {degenerated} [dump $a] ] } {
+    puts "Error: Degenerated edge is not found"
+  }
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+if { [regexp {no triangulation} [tricheck result] ] } {
+  puts "Error: Result contains not triangulated face"
+}
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b4a

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+puts "TODO OCC29883 ALL: Error :  is WRONG because number of WIRE entities in shape \"result\" is 10"
+puts "TODO OCC29883 ALL: Error :  is WRONG because number of FACE entities in shape \"result\" is 10"
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate b2 0 0 0.3
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9465.07 -s 4012.74
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b4b

@@ -0,0 +1,33 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate b2 0 0 0.3
+
+tolerance b2
+
+explode b2 v
+settolerance b2_2 1.0e-7
+
+checkshape b2
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9465.07 -s 4012.74
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b5a

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+puts "TODO OCC29860 ALL: Error :  is WRONG because number of WIRE entities in shape \"result\" is 10"
+puts "TODO OCC29860 ALL: Error :  is WRONG because number of FACE entities in shape \"result\" is 10"
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate b2 0 0 2
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9465.86 -s 4013.42
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_b5b

@@ -0,0 +1,33 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate b2 0 0 2
+
+tolerance b2
+
+explode b2 v
+settolerance b2_2 1.0e-7
+
+checkshape b2
+
+bcut result b1 b2
+checkshape result
+checkprops result -v 9465.86 -s 4013.42
+checknbshapes result -face 9 -wire 9
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}
+
+vdisplay result
+vsetdispmode 1
+vviewparams -scale 52.4803 -proj 0.285421 0.0158136 0.958272 -up -0.810772 0.537155 0.232624 -at 29.6172 13.0268 3.82234 -eye 47.7439 14.0311 64.6808
+
+checkview -screenshot -3d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_i1001

@@ -0,0 +1,26 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+explode b1 f
+explode b2 f
+smallview
+don b1_5 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00015253053837904724 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i1002

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+removeloc b2 b2
+
+explode b1 f
+explode b2 f
+smallview
+don b1_5 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00039718358540697849 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i1003

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc b2 b2
+
+explode b1 f
+explode b2 f
+smallview
+don b1_5 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 5.0314111870170835e-005 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i1004

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate b2 0 0 0.3
+
+explode b1 f
+explode b2 f
+smallview
+don b1_5 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00011289757099748416 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i1005

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate b2 0 0 2
+
+explode b1 f
+explode b2 f
+smallview
+don b1_5 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 7.7125880147734232e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i1006

@@ -0,0 +1,27 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug25542_cylinder.brep] b1
+restore [locate_data_file bug25542_cone.brep] b2
+
+explode b1 f
+explode b2 f
+
+smallview
+don b1_1 b2_1
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_1 b2_1 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.0016015772839744358 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2001

@@ -0,0 +1,27 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_5
+smallview
+don f1 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.693336906196208e-008 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2002

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+removeloc b2 b2
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_5
+smallview
+don f1 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.9119212307774807e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2003

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc b2 b2
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_5
+smallview
+don f1 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.412545178522274e-008 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2004

@@ -0,0 +1,31 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate b2 0 0 0.3
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_5
+smallview
+don f1 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.6039506888710934e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2005

@@ -0,0 +1,31 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-tool.brep] b2
+
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate b2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate b2 0 0 2
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_5
+smallview
+don f1 b2_2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 6.4774617011651419e-006 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i2006

@@ -0,0 +1,27 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug25542_cylinder.brep] b1
+restore [locate_data_file bug25542_cone.brep] b2
+
+explode b1 f
+explode b2 f
+nurbsconvert f1 b1_1
+smallview
+don f1 b2_1
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 b2_1 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 6.4791886898294872e-006 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i3001

@@ -0,0 +1,25 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+explode b1 f
+smallview
+don b1_5 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00015253053837762444 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i3002

@@ -0,0 +1,28 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+removeloc f2 f2
+
+explode b1 f
+smallview
+don b1_5 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00039718358530349535 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i3003

@@ -0,0 +1,28 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc f2 f2
+
+explode b1 f
+smallview
+don b1_5 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 5.4742962483090032e-005 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i3004

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate f2 0 0 0.3
+
+explode b1 f
+smallview
+don b1_5 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 0.00011289757087827709 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i3005

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate f2 0 0 2
+
+explode b1 f
+smallview
+don b1_5 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves b1_5 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 7.7124681583892622e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i4001

@@ -0,0 +1,26 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+explode b1 f
+nurbsconvert f1 b1_5
+smallview
+don f1 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.6933365231971514e-008 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i4002

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+removeloc f2 f2
+
+explode b1 f
+nurbsconvert f1 b1_5
+smallview
+don f1 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.9119209602049977e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i4003

@@ -0,0 +1,29 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+removeloc f2 f2
+
+explode b1 f
+nurbsconvert f1 b1_5
+smallview
+don f1 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.4125458576041928e-008 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i4004

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 -5 
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 20 
+ttranslate f2 0 0 0.3
+
+explode b1 f
+nurbsconvert f1 b1_5
+smallview
+don f1 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 2.6068687218615023e-007 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i4005

@@ -0,0 +1,30 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+restore [locate_data_file bug29807-obj.brep] b1
+restore [locate_data_file bug29807-cone.brep] f2
+
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7 -7.14142842854285 0 5
+trotate f2 +23.85857157145715500000 +12.00000000000000000000 +5.50000000000000000000 7.1414284285428495 7.0000000000000009 -0 -38
+ttranslate f2 0 0 2
+
+explode b1 f
+nurbsconvert f1 b1_5
+smallview
+don f1 f2
+fit
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} [bopcurves f1 f2 -2d] full Toler NbCurv
+
+checkreal Tolerance $Toler 6.4774618953705733e-006 0.0 0.01
+
+if {$NbCurv != 2} {
+  puts "Error: Please check NbCurves for intersector"
+}
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+smallview -2D-
+2dfit
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2d.png

tests/bugs/modalg_7/bug29807_i5001

@@ -0,0 +1,70 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+puts "TODO OCC29883 ALL: Error in res_2: T=0"
+puts "TODO OCC29883 ALL: Error: 0 vertices are expected but 2 are found"
+puts "TODO OCC29883 ALL: Error :  is WRONG because number of VERTEX entities in shape \"result\" is 3"
+
+foreach a [directory res*] {unset $a}
+
+binrestore [locate_data_file bug29807_f1.bin] f1
+binrestore [locate_data_file bug29807_f2.bin] f2
+
+mksurface s1 f1
+mksurface s2 f2
+
+trim s1 s1
+trim s2 s2
+
+intersect res s1 s2 1.0e-4
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  renamevar res res_1
+}
+
+bclearobjects
+bcleartools
+
+set ic 1
+set AllowRepeat 1
+while { $AllowRepeat != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeat 0
+  } else {
+   
+    bounds res_$ic U1 U2
+    
+    if {[dval U2-U1] < 1.0e-20} {
+      puts "Error: Wrong curve's range!"
+    }
+    
+    xdistcs res_$ic s1 U1 U2 100 2.0e-7
+    xdistcs res_$ic s2 U1 U2 100 2.0e-7
+    
+    mkedge ee res_$ic
+    baddobjects ee
+    incr ic
+  }
+}
+
+bfillds
+bbuild result
+
+smallview
+don result*
+fit
+
+# Check gaps between edges in ce
+checksection result -r 0
+checkmaxtol result -min_tol 2.0e-7
+
+checknbshapes result -edge 2 -vertex 2
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_i5002

@@ -0,0 +1,60 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+foreach a [directory res*] {unset $a}
+
+cone s1 11.4307383137554 3.49999999999979 -89.7537975119388 0 0 1 1 0 0 80.000725670142287835190342147806 9.45659107381736
+cone s2 -3.6479413426839 -11.578679656441 -89.9782110643133 0 0 1 0 1 0 5 0.250951325477062
+
+intersect res s1 s2 1.0e-4
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  renamevar res res_1
+}
+
+bclearobjects
+bcleartools
+
+set ic 1
+set AllowRepeat 1
+while { $AllowRepeat != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeat 0
+  } else {
+   
+    bounds res_$ic U1 U2
+    
+    if {[dval U2-U1] < 1.0e-20} {
+      puts "Error: Wrong curve's range!"
+    }
+    
+    xdistcs res_$ic s1 U1 U2 100 2.0e-7
+    xdistcs res_$ic s2 U1 U2 100 2.0e-7
+    
+    mkedge ee res_$ic
+    baddobjects ee
+    incr ic
+  }
+}
+
+bfillds
+bbuild result
+
+smallview
+don result*
+fit
+
+# Check gaps between edges in ce
+checksection result -r 0
+checkmaxtol result -min_tol 2.0e-7
+
+checknbshapes result -edge 3 -vertex 3
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/modalg_7/bug29807_sc01

@@ -0,0 +1,42 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+# The aim of this test is to obtain the same result on
+# different platforms (Windows, Linux, MacOS etc.)
+
+pload QAcommands
+restore [locate_data_file OCC13116_sh1.brep] b1
+restore [locate_data_file OCC13116_sh2.brep] b2
+
+explode b1 f
+explode b2 f
+
+mksurface s1 b1_3
+mksurface s2 b2_1
+
+regexp {Radius of curvature is +([-0-9.+eE]+)} [OCC29807 s1 s2 1.5704836137868956 3.0501269066007808e-007 3.9658833912538207 1.5704836137865030] full R1
+if {$R1 != -1.0} {
+  puts "Error in R1 computation"
+}
+
+regexp {Radius of curvature is +([-0-9.+eE]+)} [OCC29807 s1 s2 1.6025602743734420 -6.1366790760075673e-007 0.82433854035089271 1.5390323792163476] full R2
+if {$R2 != -1.0} {
+  puts "Error in R2 computation"
+}
+
+regexp {Radius of curvature is +([-0-9.+eE]+)} [OCC29807 s1 s2 1.6375006167098363 -8.0934069046634249e-007 0.82433889938148752 1.5040920368799497] full R3
+if {$R3 != -1.0} {
+  puts "Error in R3 computation"
+}
+
+regexp {Radius of curvature is +([-0-9.+eE]+)} [OCC29807 s1 s2 1.7042049066248015 -1.1250135685259011e-006 0.82433908359918306 1.4373877469650012] full R4
+if {$R4 != -1.0} {
+  puts "Error in R4 computation"
+}
+
+regexp {Radius of curvature is +([-0-9.+eE]+)} [OCC29807 s1 s2 1.7772619860554566 -1.3863691492588259e-006 0.82433916835888565 1.3643306675343436] full R5
+if { abs($R5-19.71766721319873) > 1.0e-7} {
+  puts "Error in R5 computation"
+}

tests/bugs/modalg_7/bug29807_svm01

@@ -0,0 +1,33 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+# Faces f2 and f6 are taken from the test case "boolean volumemaker C5 (C6)"
+
+cone con_f2 0 518.47000000000003 0 0 -1 1.1102230246251565e-016 14.999999999912038 0
+mkface f2 con_f2 0 6.2831853071795862 0 1000000
+
+cone con_f6 0 -440.74363604000001 0 0 -1 1.1102230246251565e-016 45.110284878807235 0
+mkface f6 con_f6 0 6.2831853071795862 0 1000000
+
+bsection result f2 f6
+
+smallview
+don result
+fit
+disp f2 f6
+
+checksection result -r 0
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+
+checkmaxtol result -ref 7.3189259943803184e-007
+checkprops result -l 2202.91
+checknbshapes result -vertex 1 -edge 1
+
+checkshape result
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}

tests/bugs/modalg_7/bug29807_svm02

@@ -0,0 +1,33 @@
+puts "========"
+puts "0029807: Impossible to cut cone from prism"
+puts "========"
+puts ""
+
+# Faces f2 and f6 are taken from the test case "boolean volumemaker E7"
+
+cone con_f1 0 -60.919306349999999 0 0 -1 1.1102230246251565e-016 28.800000000062262 0
+mkface f1 con_f1 0 6.2831853071795862 0 1000000
+
+cone con_f5 0 -309.47272469000001 0 0 -1 1.1102230246251565e-016 43.999999999485127 0
+mkface f5 con_f5 0 6.2831853071795862 0 1000000
+
+bsection result f1 f5
+
+smallview
+don result
+fit
+disp f1 f5
+
+checksection result -r 0
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_3d.png
+
+checkmaxtol result -ref 6.6226289034767669e-007
+checkprops result -l 1993.34
+checknbshapes result -vertex 1 -edge 1
+
+checkshape result
+
+if {[regexp "Faulties" [bopargcheck result]]} {
+  puts "Error: bopargcheck has found some faulties in result"
+}

tests/bugs/modalg_7/bug29824

@@ -0,0 +1,54 @@
+puts "========"
+puts "0029824: Intersection of cylinder and sphere is incorrect"
+puts "========"
+puts ""
+
+restore [locate_data_file OCC20964_revsolid.brep] b1 
+restore [locate_data_file OCC20964_sphere.brep] b2 
+explode b1 f
+explode b2 f
+donly b1_7 b2_1
+
+bop b1_7 b2_1
+bopsection result
+
+checknbshapes result -edge 3
+
+foreach a [explode result e] {
+  mkcurve cc $a
+  bounds cc U1 U2
+  
+  if {[dval U2-U1] < 1.0e-9} {
+    puts "Error: Wrong curve's range!"
+  }
+  
+  set aStep [expr [dval U2-U1]/100.0 ]
+  set isFirst 1
+  set aCosPrev 1.0
+  dset dx1 0.0
+  dset dy1 0.0
+  dset dz1 0.0
+  for {set aU [dval U1]} {$aU <= [dval U2]} {set aU [expr $aU + $aStep]} {
+    cvalue cc $aU xx yy zz dx dy dz
+    
+    if {!$isFirst} {
+      set m1 [module dx1 dy1 dz1]
+      set m2 [module dx dy dz]
+      set aCos [dval dx1*dx+dy1*dy+dz1*dz]
+      set aCos [expr $aCos/($m1*$m2)]
+      
+      if {abs([expr $aCos - $aCosPrev]) > 0.05} {
+        puts "Error: It seems that the curve $a change it direction at the point $aU. Please recheck."
+        break
+      }
+      
+      set aCosPrev $aCos
+    }
+    
+    set isFirst 0
+    
+    dset dx1 dx
+    dset dy1 dy
+    dset dz1 dz
+  }
+}

tests/bugs/moddata_2/bug565

@@ -1,9 +1,7 @@
-puts "========| OCC565 |========"
+puts "========"
-##################################
+puts "OCC565: Can not intersect two trimmed conical surfaces"
-## Can not intersect two trimmed conical surfaces 
+puts "========"
-##################################
+puts ""
-
-puts "TODO OCC28016 Linux: Error: 1 is expected but .* is found!"
 
 set GoodNbCurv 1