0024053: Section between plane and sphere is not correct
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
index a70db80..e546edd 100755 (executable)
@@ -1,7 +1,22 @@
-// File:      IntTools_FaceFace.cxx
-// Created:   Thu Nov 23 14:52:53 2000
-// Author:    Michael KLOKOV
-// Copyright: OPEN CASCADE 2000
+// Created on: 2000-11-23
+// Created by: Michael KLOKOV
+// Copyright (c) 2000-2012 OPEN CASCADE SAS
+//
+// The content of this file is subject to the Open CASCADE Technology Public
+// License Version 6.5 (the "License"). You may not use the content of this file
+// except in compliance with the License. Please obtain a copy of the License
+// at http://www.opencascade.org and read it completely before using this file.
+//
+// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
+// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
+//
+// The Original Code and all software distributed under the License is
+// distributed on an "AS IS" basis, without warranty of any kind, and the
+// Initial Developer hereby disclaims all such warranties, including without
+// limitation, any warranties of merchantability, fitness for a particular
+// purpose or non-infringement. Please see the License for the specific terms
+// and conditions governing the rights and limitations under the License.
+
 
 
 #include <IntTools_FaceFace.ixx>
 #include <GeomProjLib.hxx>
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
-//
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
 #include <TopExp_Explorer.hxx>
 #include <BRepTools.hxx>
 #include <BRepAdaptor_Surface.hxx>
 
-#include <BOPTColStd_Dump.hxx>
-
 #include <IntTools_Curve.hxx>
 #include <IntTools_Tools.hxx>
 #include <IntTools_Tools.hxx>
 #include <IntTools_TopolTool.hxx>
 #include <IntTools_PntOnFace.hxx>
 #include <IntTools_PntOn2Faces.hxx>
-#include <IntTools_Context.hxx>
+#include <BOPInt_Context.hxx>
 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
 
+static
+  void RefineVector(gp_Vec2d& aV2D);
 
-//modified by NIZNHY-PKV Fri Nov 25 12:03:55 2011f
 static
   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
-//modified by NIZNHY-PKV Fri Nov 25 12:03:58 2011t        
 //
 static
   void TolR3d(const TopoDS_Face& ,
@@ -178,7 +190,8 @@ static
                                        const IntTools_LineConstructor&                theLConstructor,
                                        const Standard_Boolean                         theAvoidLConstructor,
                                        IntPatch_SequenceOfLine&                       theNewLines,
-                                       Standard_Real&                                 theReachedTol3d);
+                                       Standard_Real&                                 theReachedTol3d,
+                                       const Handle(BOPInt_Context)& );
 
 static 
   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
@@ -187,7 +200,8 @@ static
                                          const TopoDS_Face&        theFace2,
                                          const Standard_Real       theOtherParameter,
                                          const Standard_Boolean    bIncreasePar,
-                                         Standard_Real&            theNewParameter);
+                                         Standard_Real&            theNewParameter,
+                                         const Handle(BOPInt_Context)& );
 
 static 
   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
@@ -215,7 +229,8 @@ static
                                       const TopoDS_Face&                   theFace2,
                                       Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
                                       Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
-                                      Handle(TColStd_HArray1OfReal)&       theResultRadius);
+                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
+                                      const Handle(BOPInt_Context)& );
 
 static
   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
@@ -236,12 +251,12 @@ static
                                   Handle(GeomAdaptor_HSurface) theGASurface);
 
 static
-gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
-                          const gp_Pnt2d&     theOriginalPoint,
-                          Handle(GeomAdaptor_HSurface) theGASurface);
+  gp_Pnt2d AdjustByNeighbour(const gp_Pnt2d&     theaNeighbourPoint,
+                            const gp_Pnt2d&     theOriginalPoint,
+                            Handle(GeomAdaptor_HSurface) theGASurface);
 static
-Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
-                                   const gp_Sphere& theSph);
+  Standard_Boolean  ApproxWithPCurves(const gp_Cylinder& theCyl, 
+                                     const gp_Sphere& theSph);
 
 static void  PerformPlanes(const Handle(GeomAdaptor_HSurface)& theS1, 
                           const Handle(GeomAdaptor_HSurface)& theS2, 
@@ -262,6 +277,7 @@ static
   void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
                        const Handle(GeomAdaptor_HSurface)& aHS2,
                        Standard_Integer& iDegMin,
+                       Standard_Integer& iNbIter,
                        Standard_Integer& iDegMax);
 
 static
@@ -279,12 +295,42 @@ static
   Standard_Integer IndexType(const GeomAbs_SurfaceType aType);
 
 //
+static
+  Standard_Real MaxSquareDistance (const Standard_Real aT,
+                                  const Handle(Geom_Curve)& aC3D,
+                                  const Handle(Geom2d_Curve)& aC2D1,
+                                  const Handle(Geom2d_Curve)& aC2D2,
+                                  const Handle(GeomAdaptor_HSurface) myHS1,
+                                  const Handle(GeomAdaptor_HSurface) myHS2,
+                                  const TopoDS_Face& aF1,
+                                  const TopoDS_Face& aF2,
+                                  const Handle(BOPInt_Context)& aCtx);
+
+static
+  Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
+                               const TopoDS_Face& aFace);
+
+//
+static
+  Standard_Real FindMaxSquareDistance (const Standard_Real aA,
+                                      const Standard_Real aB,
+                                      const Standard_Real aEps,
+                                      const Handle(Geom_Curve)& aC3D,
+                                      const Handle(Geom2d_Curve)& aC2D1,
+                                      const Handle(Geom2d_Curve)& aC2D2,
+                                      const Handle(GeomAdaptor_HSurface)& myHS1,
+                                      const Handle(GeomAdaptor_HSurface)& myHS2,
+                                      const TopoDS_Face& aF1,
+                                      const TopoDS_Face& aF2,
+                                      const Handle(BOPInt_Context)& aCtx);
+
 //=======================================================================
 //function : 
 //purpose  : 
 //=======================================================================
-  IntTools_FaceFace::IntTools_FaceFace()
+IntTools_FaceFace::IntTools_FaceFace()
 {
+  myIsDone=Standard_False;
   myTangentFaces=Standard_False;
   //
   myHS1 = new GeomAdaptor_HSurface ();
@@ -292,30 +338,45 @@ static
   myTolReached2d=0.; 
   myTolReached3d=0.;
   SetParameters(Standard_True, Standard_True, Standard_True, 1.e-07);
+  
+}
+//=======================================================================
+//function : SetContext
+//purpose  : 
+//======================================================================= 
+void IntTools_FaceFace::SetContext(const Handle(BOPInt_Context)& aContext)
+{
+  myContext=aContext;
+}
+//=======================================================================
+//function : Context
+//purpose  : 
+//======================================================================= 
+const Handle(BOPInt_Context)& IntTools_FaceFace::Context()const
+{
+  return myContext;
 }
 //=======================================================================
 //function : Face1
 //purpose  : 
 //======================================================================= 
-  const TopoDS_Face&  IntTools_FaceFace::Face1() const
+const TopoDS_Face& IntTools_FaceFace::Face1() const
 {
   return myFace1;
 }
-
 //=======================================================================
 //function : Face2
 //purpose  : 
 //======================================================================= 
-  const TopoDS_Face&  IntTools_FaceFace::Face2() const
+const TopoDS_Face& IntTools_FaceFace::Face2() const
 {
   return myFace2;
 }
-
 //=======================================================================
 //function : TangentFaces
 //purpose  : 
 //======================================================================= 
-  Standard_Boolean IntTools_FaceFace::TangentFaces() const
+Standard_Boolean IntTools_FaceFace::TangentFaces() const
 {
   return myTangentFaces;
 }
@@ -323,7 +384,7 @@ static
 //function : Points
 //purpose  : 
 //======================================================================= 
-  const  IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
+const IntTools_SequenceOfPntOn2Faces& IntTools_FaceFace::Points() const
 {
   return myPnts;
 }
@@ -331,7 +392,7 @@ static
 //function : IsDone
 //purpose  : 
 //======================================================================= 
-  Standard_Boolean IntTools_FaceFace::IsDone() const
+Standard_Boolean IntTools_FaceFace::IsDone() const
 {
   return myIsDone;
 }
@@ -339,7 +400,7 @@ static
 //function : TolReached3d
 //purpose  : 
 //=======================================================================
-  Standard_Real IntTools_FaceFace::TolReached3d() const
+Standard_Real IntTools_FaceFace::TolReached3d() const
 {
   return myTolReached3d;
 }
@@ -347,18 +408,18 @@ static
 //function : Lines
 //purpose  : return lines of intersection
 //=======================================================================
-  const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
+const IntTools_SequenceOfCurves& IntTools_FaceFace::Lines() const
 {
-  StdFail_NotDone_Raise_if(!myIsDone,
-                          "IntTools_FaceFace::Lines() => !myIntersector.IsDone()");
+  StdFail_NotDone_Raise_if
+    (!myIsDone,
+     "IntTools_FaceFace::Lines() => !myIntersector.IsDone()");
   return mySeqOfCurve;
 }
-
 //=======================================================================
 //function : TolReached2d
 //purpose  : 
 //=======================================================================
-  Standard_Real IntTools_FaceFace::TolReached2d() const
+Standard_Real IntTools_FaceFace::TolReached2d() const
 {
   return myTolReached2d;
 }
@@ -366,10 +427,10 @@ static
 // function: SetParameters
 //
 // =======================================================================
-  void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
-                                       const Standard_Boolean ToApproxC2dOnS1,
-                                       const Standard_Boolean ToApproxC2dOnS2,
-                                       const Standard_Real ApproximationTolerance) 
+void IntTools_FaceFace::SetParameters(const Standard_Boolean ToApproxC3d,
+                                     const Standard_Boolean ToApproxC2dOnS1,
+                                     const Standard_Boolean ToApproxC2dOnS2,
+                                     const Standard_Real ApproximationTolerance) 
 {
   myApprox = ToApproxC3d;
   myApprox1 = ToApproxC2dOnS1;
@@ -380,12 +441,10 @@ static
 //function : SetList
 //purpose  : 
 //=======================================================================
-
 void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
 {
   myListOfPnts = aListOfPnts;  
 }
-
 //=======================================================================
 //function : Perform
 //purpose  : intersect surfaces of the faces
@@ -403,6 +462,10 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
   Handle(IntTools_TopolTool) dom1, dom2;
   BRepAdaptor_Surface aBAS1, aBAS2;
   //
+  if (myContext.IsNull()) {
+    myContext=new BOPInt_Context;
+  }
+  //
   mySeqOfCurve.Clear();
   myTolReached2d=0.;
   myTolReached3d=0.;
@@ -462,38 +525,35 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                  mySeqOfCurve, myTangentFaces);
 
     myIsDone = Standard_True;
-
-    if(myTangentFaces) {
-      return;
-    }
-    //
-    NbLinPP = mySeqOfCurve.Length();
-    if(NbLinPP == 0) {
-      return;
-    }
-
-    Standard_Real aTolFMax;
-    //
-    myTolReached3d = 1.e-7;
-    //
-    aTolFMax=Max(aTolF1, aTolF2);
-    //
-    if (aTolFMax>myTolReached3d) {
-      myTolReached3d=aTolFMax;
-    }
-    myTolReached2d = myTolReached3d;
-    //
-    if (bReverse) {
-      Handle(Geom2d_Curve) aC2D1, aC2D2;
+    
+    if(!myTangentFaces) {
       //
-      aNbLin=mySeqOfCurve.Length();
-      for (i=1; i<=aNbLin; ++i) {
-       IntTools_Curve& aIC=mySeqOfCurve(i);
-       aC2D1=aIC.FirstCurve2d();
-       aC2D2=aIC.SecondCurve2d();
+      NbLinPP = mySeqOfCurve.Length();
+      if(NbLinPP) {
+       Standard_Real aTolFMax;
        //
-       aIC.SetFirstCurve2d(aC2D2);
-       aIC.SetSecondCurve2d(aC2D1);
+       myTolReached3d = 1.e-7;
+       //
+       aTolFMax=Max(aTolF1, aTolF2);
+       //
+       if (aTolFMax>myTolReached3d) {
+         myTolReached3d=aTolFMax;
+       }
+       myTolReached2d = myTolReached3d;
+       //
+       if (bReverse) {
+         Handle(Geom2d_Curve) aC2D1, aC2D2;
+         //
+         aNbLin=mySeqOfCurve.Length();
+         for (i=1; i<=aNbLin; ++i) {
+           IntTools_Curve& aIC=mySeqOfCurve(i);
+           aC2D1=aIC.FirstCurve2d();
+           aC2D2=aIC.SecondCurve2d();
+           //
+           aIC.SetFirstCurve2d(aC2D2);
+           aIC.SetSecondCurve2d(aC2D1);
+         }
+       }
       }
     }
     return;
@@ -688,39 +748,85 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
   GeomAbs_SurfaceType aType1, aType2;
   //
   aNbLin=myIntersector.NbLines();
+  if (!aNbLin) {
+    return;
+  }
   //
   aType1=myHS1->Surface().GetType();
   aType2=myHS2->Surface().GetType();
   //
-  if (aNbLin==2 &&
-      aType1==GeomAbs_Cylinder &&
-      aType2==GeomAbs_Cylinder) {
-    Handle(IntPatch_Line) aIL1, aIL2;
-    IntPatch_IType aTL1, aTL2;
-    //
-    aIL1=myIntersector.Line(1);
-    aIL2=myIntersector.Line(2);
-    aTL1=aIL1->ArcType();
-    aTL2=aIL2->ArcType();
-    if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
-      Standard_Real aD, aDTresh, dTol;
-      gp_Lin aL1, aL2;
-      //
-      dTol=1.e-8;
-      aDTresh=1.5e-6;
+  if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
+    if (aNbLin==2){ 
+      Handle(IntPatch_Line) aIL1, aIL2;
+      IntPatch_IType aTL1, aTL2;
       //
-      aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
-      aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
-      aD=aL1.Distance(aL2);
-      aD=0.5*aD;
-      if (aD<aDTresh) {
-       myTolReached3d=aD+dTol;
+      aIL1=myIntersector.Line(1);
+      aIL2=myIntersector.Line(2);
+      aTL1=aIL1->ArcType();
+      aTL2=aIL2->ArcType();
+      if (aTL1==IntPatch_Lin && aTL2==IntPatch_Lin) {
+       Standard_Real aD, aDTresh, dTol;
+       gp_Lin aL1, aL2;
+       //
+       dTol=1.e-8;
+       aDTresh=1.5e-6;
+       //
+       aL1=Handle(IntPatch_GLine)::DownCast(aIL1)->Line();
+       aL2=Handle(IntPatch_GLine)::DownCast(aIL2)->Line();
+       aD=aL1.Distance(aL2);
+       aD=0.5*aD;
+       if (aD<aDTresh) {
+         myTolReached3d=aD+dTol;
+       }
+       return;
       }
     }
-  }
+    //ZZ
+    if (aNbLin) {// Check the distances
+      Standard_Integer i, aNbP, j ;
+      Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
+      //
+      aD2Max=0.;
+      aNbP=10;
+      aNbLin=mySeqOfCurve.Length();
+      //
+      for (i=1; i<=aNbLin; ++i) {
+       const IntTools_Curve& aIC=mySeqOfCurve(i);
+       const Handle(Geom_Curve)& aC3D=aIC.Curve();
+       const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
+       const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
+       //
+       if (aC3D.IsNull()) {
+         continue;
+       }
+       const Handle(Geom_BSplineCurve)& aBC=
+         Handle(Geom_BSplineCurve)::DownCast(aC3D);
+       if (aBC.IsNull()) {
+         continue;
+       }
+       //
+       aT1=aBC->FirstParameter();
+       aT2=aBC->LastParameter();
+       //
+       aEps=0.01*(aT2-aT1);
+       dT=(aT2-aT1)/aNbP;
+       for (j=1; j<aNbP; ++j) {
+         aT11=aT1+j*dT;
+         aT12=aT11+dT;
+         aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
+                                   myHS1, myHS2, myFace1, myFace2, myContext);
+         if (aD2>aD2Max) {
+           aD2Max=aD2;
+         }
+       }
+      }//for (i=1; i<=aNbLin; ++i) {
+      //
+      myTolReached3d=sqrt(aD2Max);
+    }// if (aNbLin) 
+  }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
+  //
   //904/G3 f
-  if (aType1==GeomAbs_Plane &&
-      aType2==GeomAbs_Plane) {
+  else if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
     Standard_Real aTolF1, aTolF2, aTolFMax, aTolTresh;
     //
     aTolTresh=1.e-7;
@@ -732,11 +838,11 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
     if (aTolFMax>aTolTresh) {
       myTolReached3d=aTolFMax;
     }
-  }
+  }//if (aType1==GeomAbs_Plane && aType2==GeomAbs_Plane) {
   //t
   //IFV Bug OCC20297 
-  if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
-     (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
+  else if((aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ||
+         (aType2 == GeomAbs_Cylinder && aType1 == GeomAbs_Plane)) {
     if(aNbLin == 1) {
       const Handle(IntPatch_Line)& aIL1 = myIntersector.Line(1);
       if(aIL1->ArcType() == IntPatch_Circle) {
@@ -756,8 +862,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
        Standard_Real eps = 1.e-14;
        if(!aPlDir.IsEqual(aCirDir, eps)) {
          Standard_Integer aNbP = 11;
-         Standard_Real dt = 2.*PI / (aNbP - 1), t;
-         for(t = 0.; t < 2.*PI; t += dt) {
+         Standard_Real dt = 2.*M_PI / (aNbP - 1), t;
+         for(t = 0.; t < 2.*M_PI; t += dt) {
            Standard_Real d = aPln.Distance(ElCLib::Value(t, aCir)); 
            if(myTolReached3d < d) myTolReached3d = d;
          }
@@ -765,11 +871,11 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
        }
       } //aIL1->ArcType() == IntPatch_Circle
     } //aNbLin == 1
-  } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) ...
+  } // aType1 == GeomAbs_Cylinder && aType2 == GeomAbs_Plane) 
   //End IFV Bug OCC20297
   //
-  if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
-      (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
+  else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
+          (aType2==GeomAbs_Plane && aType1==GeomAbs_Torus)) {
     aNbLin=mySeqOfCurve.Length();
     if (aNbLin!=1) {
       return;
@@ -784,7 +890,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
     //
     const IntTools_Curve& aIC=mySeqOfCurve(1);
     const Handle(Geom_Curve)& aC3D=aIC.Curve();
-    const Handle(Geom_BSplineCurve)& aBS=Handle(Geom_BSplineCurve)::DownCast(aC3D);
+    const Handle(Geom_BSplineCurve)& aBS=
+      Handle(Geom_BSplineCurve)::DownCast(aC3D);
     if (aBS.IsNull()) {
       return;
     }
@@ -827,18 +934,13 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
     }
   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
   //
-  if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
-      (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
-    Standard_Boolean bIsDone;
+  else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
+          (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
     Standard_Integer i, j, aNbP;
-    Standard_Real aT, aT1, aT2, dT, aU1, aV1, aU2, aV2;
-    Standard_Real aDSmax, aDS1, aDS2, aDS;
-    gp_Pnt2d aP2D1, aP2D2;
-    gp_Pnt aP3D, aP3D1, aP3D2;
-    IntTools_Context aCtx;
+    Standard_Real aT, aT1, aT2, dT, aD2max, aD2;
     //
     aNbLin=mySeqOfCurve.Length();
-    aDSmax=-1.;
+    aD2max=0.;
     aNbP=11;
     //
     for (i=1; i<=aNbLin; ++i) {
@@ -850,7 +952,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
       if (aC3D.IsNull()) {
        continue;
       }
-      const Handle(Geom_BSplineCurve)& aBC=Handle(Geom_BSplineCurve)::DownCast(aC3D);
+      const Handle(Geom_BSplineCurve)& aBC=
+       Handle(Geom_BSplineCurve)::DownCast(aC3D);
       if (aBC.IsNull()) {
        return;
       }
@@ -865,59 +968,117 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
          aT=aT2;
        }
        //
-       aC3D->D0(aT, aP3D);
-       // 1
-       if (!aC2D1.IsNull()) {
-         aC2D1->D0(aT, aP2D1);
-         aP2D1.Coord(aU1, aV1);
-         myHS1->D0(aU1, aV1, aP3D1);
-         aDS1=aP3D.SquareDistance(aP3D1);
-         if (aDS1>aDSmax) {
-           aDSmax=aDS1;
-         }
-       }
-       // 2
-       if (!aC2D2.IsNull()) {
-         aC2D2->D0(aT, aP2D2);
-         aP2D2.Coord(aU2, aV2);
-         myHS2->D0(aU2, aV2, aP3D2);
-         aDS2=aP3D.SquareDistance(aP3D2);
-         if (aDS2>aDSmax) {
-           aDSmax=aDS2;
-         }
-       }
-       // 3
-       GeomAPI_ProjectPointOnSurf& aPPS1=aCtx.ProjPS(myFace1);
-       aPPS1.Perform(aP3D);
-       bIsDone=aPPS1.IsDone();
-       if (bIsDone) {
-         aPPS1.LowerDistanceParameters(aU1, aV1);
-         myHS1->D0(aU1, aV1, aP3D1);
-         aDS1=aP3D.SquareDistance(aP3D1);
-         if (aDS1>aDSmax) {
-           aDSmax=aDS1;
-         }
-       }
-       // 4
-       GeomAPI_ProjectPointOnSurf& aPPS2=aCtx.ProjPS(myFace2);
-       aPPS2.Perform(aP3D);
-       bIsDone=aPPS2.IsDone();
-       if (bIsDone) {
-         aPPS2.LowerDistanceParameters(aU2, aV2);
-         myHS2->D0(aU2, aV2, aP3D2);
-         aDS2=aP3D.SquareDistance(aP3D2);
-         if (aDS2>aDSmax) {
-           aDSmax=aDS2;
-         }
+       aD2=MaxSquareDistance(aT, aC3D, aC2D1, aC2D2,
+                             myHS1, myHS2, myFace1, myFace2, myContext);
+       if (aD2>aD2max) {
+         aD2max=aD2;
        }
       }//for (j=0; j<aNbP; ++j) {
+     
     }//for (i=1; i<=aNbLin; ++i) {
     //
-    aDS=myTolReached3d*myTolReached3d;
-    if (aDSmax > aDS) {
-      myTolReached3d=sqrt(aDSmax);
+    aD2=myTolReached3d*myTolReached3d;
+    if (aD2max > aD2) {
+      myTolReached3d=sqrt(aD2max);
     }
   }//if((aType1==GeomAbs_SurfaceOfRevolution ...
+  else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere)) {
+    Standard_Integer i, j, aNbP;
+    Standard_Real aT1, aT2, dT, aD2max, aD2, aEps, aT11, aT12;
+    //
+    aNbLin=mySeqOfCurve.Length();
+    aD2max=0.;
+    aNbP=10;
+    //
+    for (i=1; i<=aNbLin; ++i) {
+      const IntTools_Curve& aIC=mySeqOfCurve(i);
+      const Handle(Geom_Curve)& aC3D=aIC.Curve();
+      const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
+      const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
+      //
+      const Handle(Geom2d_BSplineCurve)& aBC2D1=
+        Handle(Geom2d_BSplineCurve)::DownCast(aC2D1);
+      const Handle(Geom2d_BSplineCurve)& aBC2D2=
+        Handle(Geom2d_BSplineCurve)::DownCast(aC2D2);
+      //
+      if (aBC2D1.IsNull() && aBC2D2.IsNull()) {
+        return;
+      }
+      //
+      if (!aBC2D1.IsNull()) {
+        aT1=aBC2D1->FirstParameter();
+        aT2=aBC2D1->LastParameter();
+      }
+      else {
+        aT1=aBC2D2->FirstParameter();
+        aT2=aBC2D2->LastParameter();
+      }
+      //
+      aEps=0.01*(aT2-aT1);
+      dT=(aT2-aT1)/aNbP;
+      for (j=0; j<aNbP; ++j) {
+        aT11=aT1+j*dT;
+        aT12=aT11+dT;
+        if (j==aNbP-1) {
+          aT12=aT2;
+        }
+        //
+        aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
+                                  myHS1, myHS2, myFace1, myFace2, myContext);
+        if (aD2>aD2max) {
+          aD2max=aD2;
+        }
+      }//for (j=0; j<aNbP; ++j) {
+     
+    }//for (i=1; i<=aNbLin; ++i) {
+    //
+    aD2=myTolReached3d*myTolReached3d;
+    if (aD2max > aD2) {
+      myTolReached3d=sqrt(aD2max);
+    }
+  }//else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ...
+  else if (!myApprox) {
+    Standard_Integer i, aNbP, j ;
+    Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
+    //
+    aD2Max=0.;
+    aNbLin=mySeqOfCurve.Length();
+    //
+    for (i=1; i<=aNbLin; ++i) {
+      const IntTools_Curve& aIC=mySeqOfCurve(i);
+      const Handle(Geom_Curve)& aC3D=aIC.Curve();
+      const Handle(Geom2d_Curve)& aC2D1=aIC.FirstCurve2d();
+      const Handle(Geom2d_Curve)& aC2D2=aIC.SecondCurve2d();
+      //
+      if (aC3D.IsNull()) {
+       continue;
+      }
+      const Handle(Geom_BSplineCurve)& aBC=
+       Handle(Geom_BSplineCurve)::DownCast(aC3D);
+      if (aBC.IsNull()) {
+       continue;
+      }
+      //
+      aT1=aBC->FirstParameter();
+      aT2=aBC->LastParameter();
+      //
+      aEps=0.0001*(aT2-aT1);
+      aNbP=11;
+      dT=(aT2-aT1)/aNbP;
+      for (j=1; j<aNbP-1; ++j) {
+       aT11=aT1+j*dT;
+       aT12=aT11+dT;
+       aD2=FindMaxSquareDistance(aT11, aT12, aEps, aC3D, aC2D1, aC2D2,
+                                 myHS1, myHS2, myFace1, myFace2, myContext);
+       if (aD2>aD2Max) {
+         aD2Max=aD2;
+       }
+      }
+    }//for (i=1; i<=aNbLin; ++i) {
+    myTolReached3d=sqrt(aD2Max);
+  }
+  //modified by NIZNHY-PKV Thu Aug 30 13:31:12 2012t
 }
 //=======================================================================
 //function : MakeCurve
@@ -927,7 +1088,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                                    const Handle(Adaptor3d_TopolTool)& dom1,
                                    const Handle(Adaptor3d_TopolTool)& dom2) 
 {
-  Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor;
+  Standard_Boolean bDone, rejectSurface, reApprox, bAvoidLineConstructor,
+                   bPCurvesOk;
   Standard_Boolean ok;
   Standard_Integer i, j, aNbParts;
   Standard_Real fprm, lprm;
@@ -941,6 +1103,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
   //
   rejectSurface = Standard_False;
   reApprox = Standard_False;
+  //
+  bPCurvesOk = Standard_True;
  
  reapprox:;
   
@@ -954,12 +1118,22 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
     //
     const Handle(IntPatch_WLine)& aWLine=
       Handle(IntPatch_WLine)::DownCast(L);
-    //
+    //DEBf
+    //DumpWLine(aWLine);
+    //DEBt
     anewL = ComputePurgedWLine(aWLine);
     if(anewL.IsNull()) {
       return;
     }
     L = anewL;
+    //DEBf
+    /*
+    { const Handle(IntPatch_WLine)& aWLineX=
+       Handle(IntPatch_WLine)::DownCast(L);
+      DumpWLine(aWLineX);
+    }
+    */
+    //DEBt 
     //
     if(!myListOfPnts.IsEmpty()) {
       bAvoidLineConstructor = Standard_True;
@@ -1134,7 +1308,7 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
     
     aNul=0.;
-    aPeriod=PI+PI;
+    aPeriod=M_PI+M_PI;
 
     for (i=1; i<=aNbParts; i++) {
       myLConstruct.Part(i, fprm, lprm);
@@ -1157,7 +1331,7 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
          if(P1.Distance(P2) > aTolDist) {
            Standard_Real anewpar = fprm;
 
-           if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, lprm, Standard_False, anewpar)) {
+           if(ParameterOutOfBoundary(fprm, newc, myFace1, myFace2, lprm, Standard_False, anewpar, myContext)) {
              fprm = anewpar;
            }
            aSeqFprm.Append(fprm);
@@ -1179,7 +1353,7 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
          if(P1.Distance(P2) > aTolDist) {
            Standard_Real anewpar = lprm;
 
-           if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, fprm, Standard_True, anewpar)) {
+           if(ParameterOutOfBoundary(lprm, newc, myFace1, myFace2, fprm, Standard_True, anewpar, myContext)) {
              lprm = anewpar;
            }
            aSeqFprm.Append(aNul);
@@ -1201,7 +1375,7 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
       lprm=aSeqLprm(i);
       //
       Standard_Real aRealEpsilon=RealEpsilon();
-      if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*PI) > aRealEpsilon) {
+      if (Abs(fprm) > aRealEpsilon || Abs(lprm-2.*M_PI) > aRealEpsilon) {
        //==============================================
        ////
        IntTools_Curve aCurve;
@@ -1248,14 +1422,14 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
        }
        mySeqOfCurve.Append(aCurve);
          //==============================================      
-      } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*PI) > RealEpsilon())
+      } //if (Abs(fprm) > RealEpsilon() || Abs(lprm-2.*M_PI) > RealEpsilon())
 
       else {
        //  on regarde si on garde
        //
        if (aNbParts==1) {
-//       if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*PI) < RealEpsilon()) {
-         if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*PI) <= aRealEpsilon) {
+//       if (Abs(fprm) < RealEpsilon() &&  Abs(lprm-2.*M_PI) < RealEpsilon()) {
+         if (Abs(fprm) <= aRealEpsilon && Abs(lprm-2.*M_PI) <= aRealEpsilon) {
            IntTools_Curve aCurve;
            Handle(Geom_TrimmedCurve) aTC3D=new Geom_TrimmedCurve(newc,fprm,lprm);
            aCurve.SetCurve(aTC3D);
@@ -1296,7 +1470,7 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
        //
        Standard_Real aTwoPIdiv17, u1, v1, u2, v2, Tol;
 
-       aTwoPIdiv17=2.*PI/17.;
+       aTwoPIdiv17=2.*M_PI/17.;
 
        for (j=0; j<=17; j++) {
          gp_Pnt ptref (newc->Value (j*aTwoPIdiv17));
@@ -1617,10 +1791,11 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
        if(reApprox && !rejectSurface)
          theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
        else {
-         Standard_Integer iDegMax, iDegMin;
+         Standard_Integer iDegMax, iDegMin, iNbIter;
+         //
+         ApproxParameters(myHS1, myHS2, iDegMin, iDegMax, iNbIter);
+         theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, iNbIter, Standard_True, aParType);
          //
-         ApproxParameters(myHS1, myHS2, iDegMin, iDegMax);
-         theapp3d.SetParameters(myTolApprox, tol2d, iDegMin, iDegMax, 0, Standard_True, aParType);
        }
       }
       //
@@ -1633,7 +1808,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                                           myLConstruct, 
                                           bAvoidLineConstructor, 
                                           aSeqOfL, 
-                                          aReachedTol);
+                                          aReachedTol,
+                                          myContext);
       if ( bIsDecomposited && ( myTolReached3d < aReachedTol ) )
        myTolReached3d = aReachedTol;
 
@@ -1875,7 +2051,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                    goto reapprox;
                  }
                }
-               // ###########################################
+                // ###########################################
+                bPCurvesOk = CheckPCurve(BS1, myFace2);
                aCurve.SetSecondCurve2d(BS1);
              }
              else {
@@ -1900,7 +2077,8 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                    goto reapprox;
                  }
                }
-               // ###########################################
+                // ###########################################
+                bPCurvesOk = bPCurvesOk && CheckPCurve(BS2, myFace1);
                aCurve.SetFirstCurve2d(BS2);
              }
              else { 
@@ -1909,7 +2087,35 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
                aCurve.SetFirstCurve2d(H1);
              }
              //
-             mySeqOfCurve.Append(aCurve);
+              //if points of the pcurves are out of the faces bounds
+              //create 3d and 2d curves without approximation
+              if (!bPCurvesOk) {
+                Handle(Geom2d_BSplineCurve) H1, H2;
+               bPCurvesOk = Standard_True;
+                //       
+                Handle(Geom_Curve) aBSp=MakeBSpline(WL,ifprm, ilprm);
+                
+                if(myApprox1) {
+                  H1 = MakeBSpline2d(WL, ifprm, ilprm, Standard_True);
+                 bPCurvesOk = CheckPCurve(H1, myFace1);
+                }
+                
+                if(myApprox2) {
+                  H2 = MakeBSpline2d(WL, ifprm, ilprm, Standard_False);
+                 bPCurvesOk = bPCurvesOk && CheckPCurve(H2, myFace2);
+                }
+                //
+               //if pcurves created without approximation are out of the 
+               //faces bounds, use approximated 3d and 2d curves
+               if (bPCurvesOk) {
+                 IntTools_Curve aIC(aBSp, H1, H2);
+                 mySeqOfCurve.Append(aIC);
+               } else {
+                 mySeqOfCurve.Append(aCurve);
+               }
+              } else {
+                mySeqOfCurve.Append(aCurve);
+              }
            }
            else { 
              const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
@@ -2107,9 +2313,6 @@ void IntTools_FaceFace::SetList(IntSurf_ListOfPntOn2S& aListOfPnts)
       }
     }
   }
-  if (C2d.IsNull()) {
-    BOPTColStd_Dump::PrintMessage("BuildPCurves()=> Echec ProjLib\n");
-  }
 }
 
 //=======================================================================
@@ -2348,6 +2551,15 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
     enlarge=Standard_True;
   }
   //
+  if (aType==GeomAbs_Sphere) {
+    Standard_Real dV;
+    //
+    dV=thevmax-thevmin;
+    if (dV+delta<M_PI) {
+      enlarge=Standard_True;
+    }
+  }
+  //
   if(!isuperiodic && enlarge) {
 
     if((theumin - uinf) > delta )
@@ -2507,12 +2719,12 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
 // The block is dedicated to determine whether WLine [ifprm, ilprm]
 // crosses the degenerated zone on each given surface or not.
 // If Yes -> We will not use info about surfaces during approximation
-// because inside degenerated zone of the surface the approx. alogo.
+// because inside degenerated zone of the surface the approx. algo.
 // uses wrong values of normal, etc., and resulting curve will have
 // oscillations that we would not like to have. 
-//                                               PKV Tue Feb 12 2002  
  
 
+
 static
   Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
                                     const Handle(Geom_Surface)& aS,
@@ -2709,16 +2921,23 @@ Standard_Boolean IsDegeneratedZone(const gp_Pnt2d& aP2d,
 //           of the points is less than 2.
 //=========================================================================
 Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine) {
+  Standard_Integer i, k, v, nb, nbvtx;
   Handle(IntPatch_WLine) aResult;
+  nbvtx = theWLine->NbVertex();
+  nb = theWLine->NbPnts();
+  if (nb==2) {
+    const IntSurf_PntOn2S& p1 = theWLine->Point(1);
+    const IntSurf_PntOn2S& p2 = theWLine->Point(2);
+    if(p1.Value().IsEqual(p2.Value(), gp::Resolution())) {
+      return aResult;
+    }
+  }
+  //
   Handle(IntPatch_WLine) aLocalWLine;
   Handle(IntPatch_WLine) aTmpWLine = theWLine;
-
   Handle(IntSurf_LineOn2S) aLineOn2S = new IntSurf_LineOn2S();
   aLocalWLine = new IntPatch_WLine(aLineOn2S, Standard_False);
-  Standard_Integer i, k, v, nb, nbvtx;
-  nbvtx = theWLine->NbVertex();
-  nb = theWLine->NbPnts();
-
   for(i = 1; i <= nb; i++) {
     aLineOn2S->Add(theWLine->Point(i));
   }
@@ -2733,10 +2952,7 @@ Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine
     nb = aLineOn2S->NbPoints();
     anEndIndex = (anEndIndex > nb) ? nb : anEndIndex;
 
-    //modified by NIZNHY-PKV Fri Nov 25 13:02:40 2011f
     if((aStartIndex > nb) || (anEndIndex <= 1)) {
-    //if((aStartIndex >= nb) || (anEndIndex <= 1)) {
-    //modified by NIZNHY-PKV Fri Nov 25 13:02:47 2011t
       continue;
     }
     k = aStartIndex;
@@ -2892,6 +3108,7 @@ Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
     }
     gp_Vec2d anormvec = aVec;
     anormvec.Normalize();
+    RefineVector(anormvec);
     Standard_Real adot1 = anormvec.Dot(anOtherVecNormal);
 
     if(fabs(adot1) < Precision::Angular())
@@ -3036,8 +3253,6 @@ Standard_Boolean CheckTangentZonesExist( const Handle(GeomAdaptor_HSurface)& the
       ( theSurface2->GetType() != GeomAbs_Torus ) )
     return Standard_False;
 
-  IntTools_Context aContext;
-
   gp_Torus aTor1 = theSurface1->Torus();
   gp_Torus aTor2 = theSurface2->Torus();
 
@@ -3064,12 +3279,13 @@ Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSur
                                     const TopoDS_Face&                   theFace2,
                                     Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
                                     Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
-                                    Handle(TColStd_HArray1OfReal)&       theResultRadius) {
+                                    Handle(TColStd_HArray1OfReal)&       theResultRadius,
+                                    const Handle(BOPInt_Context)& aContext)
+{
   Standard_Integer aResult = 0;
   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
     return aResult;
 
-  IntTools_Context aContext;
 
   TColgp_SequenceOfPnt2d aSeqResultS1, aSeqResultS2;
   TColStd_SequenceOfReal aSeqResultRad;
@@ -3125,7 +3341,7 @@ Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSur
 
     GeomAdaptor_Curve aC1( new Geom_Circle(aCircle1) );
     GeomAdaptor_Curve aC2( new Geom_Circle(aCircle2) );
-    Extrema_ExtCC anExtrema(aC1, aC2, 0, 2.*Standard_PI, 0, 2.*Standard_PI, 
+    Extrema_ExtCC anExtrema(aC1, aC2, 0, 2. * M_PI, 0, 2. * M_PI, 
                            Precision::PConfusion(), Precision::PConfusion());
        
     if ( anExtrema.IsDone() ) {
@@ -3143,7 +3359,8 @@ Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSur
 
        Standard_Integer surfit = 0;
        for ( surfit = 0; surfit < 2; surfit++ ) {
-         GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aContext.ProjPS(theFace1) : aContext.ProjPS(theFace2);
+         GeomAPI_ProjectPointOnSurf& aProjector = 
+           (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
 
          gp_Pnt aP3d = (surfit == 0) ? P1.Value() : P2.Value();
          aProjector.Perform(aP3d);
@@ -3170,7 +3387,7 @@ Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSur
          aSeqResultRad.Append( aCriteria );
 
          // torus is u and v periodic
-         const Standard_Real twoPI = Standard_PI + Standard_PI;
+         const Standard_Real twoPI = M_PI + M_PI;
          Standard_Real arr1tmp[2] = {pr1.X(), pr1.Y()};
          Standard_Real arr2tmp[2] = {pr2.X(), pr2.Y()};
 
@@ -3274,7 +3491,9 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
                                      const IntTools_LineConstructor&                theLConstructor,
                                      const Standard_Boolean                         theAvoidLConstructor,
                                      IntPatch_SequenceOfLine&                       theNewLines,
-                                     Standard_Real&                                 theReachedTol3d) {
+                                     Standard_Real&                                 theReachedTol3d,
+                                     const Handle(BOPInt_Context)& aContext) 
+{
 
   Standard_Boolean bRet, bAvoidLineConstructor;
   Standard_Integer aNbPnts, aNbParts;
@@ -3299,13 +3518,12 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
   TColStd_Array1OfListOfInteger anArrayOfLines(1, aNbPnts); 
   TColStd_Array1OfInteger       anArrayOfLineType(1, aNbPnts);
   TColStd_ListOfInteger aListOfPointIndex;
-  IntTools_Context aContext;
   
   Handle(TColgp_HArray1OfPnt2d) aTanZoneS1;
   Handle(TColgp_HArray1OfPnt2d) aTanZoneS2;
   Handle(TColStd_HArray1OfReal) aTanZoneRadius;
   Standard_Integer aNbZone = ComputeTangentZones( theSurface1, theSurface2, theFace1, theFace2,
-                                                aTanZoneS1, aTanZoneS2, aTanZoneRadius );
+                                                aTanZoneS1, aTanZoneS2, aTanZoneRadius, aContext);
   
   //
   nblines=0;
@@ -3622,7 +3840,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
                else {
                  Standard_Real anAngle = aNewVec.Angle(aVecOld);
 
-                 if((fabs(anAngle) < (Standard_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
+                 if((fabs(anAngle) < (M_PI * 0.25)) && (aNewVec.Dot(aVecOld) > 0.)) {
 
                    if(bCheckAngle1) {
                      Standard_Real U1, U2, V1, V2;
@@ -3743,7 +3961,8 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
          if(found) {
            // check point
            Standard_Real aCriteria = BRep_Tool::Tolerance(theFace1) + BRep_Tool::Tolerance(theFace2);
-           GeomAPI_ProjectPointOnSurf& aProjector = (surfit == 0) ? aContext.ProjPS(theFace2) : aContext.ProjPS(theFace1);
+           GeomAPI_ProjectPointOnSurf& aProjector = 
+             (surfit == 0) ? aContext->ProjPS(theFace2) : aContext->ProjPS(theFace1);
            Handle(GeomAdaptor_HSurface) aSurface = (surfit == 0) ? theSurface1 : theSurface2;
 
            Handle(GeomAdaptor_HSurface) aSurfaceOther = (surfit == 0) ? theSurface2 : theSurface1;
@@ -3779,7 +3998,8 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
                  foundV = ( foundV < vmin ) ? vmin : foundV;
                  foundV = ( foundV > vmax ) ? vmax : foundV;
 
-                 GeomAPI_ProjectPointOnSurf& aProjector2 = (surfit == 0) ? aContext.ProjPS(theFace1) : aContext.ProjPS(theFace2);
+                 GeomAPI_ProjectPointOnSurf& aProjector2 = 
+                   (surfit == 0) ? aContext->ProjPS(theFace1) : aContext->ProjPS(theFace2);
 
                  aP3d = aSurfaceOther->Value(foundU, foundV);
                  aProjector2.Perform(aP3d);
@@ -3996,11 +4216,12 @@ Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
                                        const TopoDS_Face&        theFace2,
                                        const Standard_Real       theOtherParameter,
                                        const Standard_Boolean    bIncreasePar,
-                                       Standard_Real&            theNewParameter) {
+                                       Standard_Real&            theNewParameter,
+                                       const Handle(BOPInt_Context)& aContext)
+{
   Standard_Boolean bIsComputed = Standard_False;
   theNewParameter = theParameter;
 
-  IntTools_Context aContext;
   Standard_Real acurpar = theParameter;
   TopAbs_State aState = TopAbs_ON;
   Standard_Integer iter = 0;
@@ -4031,7 +4252,7 @@ Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
 
     if(aPrj1.IsDone()) {
       aPrj1.LowerDistanceParameters(U, V);
-      aState = aContext.StatePointFace(theFace1, gp_Pnt2d(U, V));
+      aState = aContext->StatePointFace(theFace1, gp_Pnt2d(U, V));
     }
 
     if(aState != TopAbs_ON) {
@@ -4039,7 +4260,7 @@ Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
                
       if(aPrj2.IsDone()) {
        aPrj2.LowerDistanceParameters(U, V);
-       aState = aContext.StatePointFace(theFace2, gp_Pnt2d(U, V));
+       aState = aContext->StatePointFace(theFace2, gp_Pnt2d(U, V));
       }
     }
 
@@ -4388,11 +4609,14 @@ Standard_Boolean ClassifyLin2d(const Handle(GeomAdaptor_HSurface)& theS,
 void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
                      const Handle(GeomAdaptor_HSurface)& aHS2,
                      Standard_Integer& iDegMin,
-                     Standard_Integer& iDegMax)
+                     Standard_Integer& iDegMax,
+                     Standard_Integer& iNbIter)
+
 {
   GeomAbs_SurfaceType aTS1, aTS2;
   
   //
+  iNbIter=0;
   iDegMin=4;
   iDegMax=8;
   //
@@ -4422,6 +4646,9 @@ void ApproxParameters(const Handle(GeomAdaptor_HSurface)& aHS1,
       iDegMax=6;
     }
   }
+  if (aTS1==GeomAbs_Cylinder && aTS2==GeomAbs_Cylinder) {
+    iNbIter=1; 
+  }
 }
 //=======================================================================
 //function : Tolerances
@@ -4527,8 +4754,6 @@ Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
   } 
   return aIndex;
 }
-
-//modified by NIZNHY-PKV Fri Nov 25 11:58:12 2011f
 //=======================================================================
 //function : DumpWLine
 //purpose  : 
@@ -4538,6 +4763,7 @@ void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
   Standard_Integer i, aNbPnts; 
   Standard_Real aX, aY, aZ, aU1, aV1, aU2, aV2;
   //
+  printf(" *WLine\n");
   aNbPnts=aWLine->NbPnts();
   for (i=1; i<=aNbPnts; ++i) {
     const IntSurf_PntOn2S aPntOn2S=aWLine->Point(i);
@@ -4545,9 +4771,208 @@ void DumpWLine(const Handle(IntPatch_WLine)& aWLine)
     aP3D.Coord(aX, aY, aZ);
     aPntOn2S.Parameters(aU1, aV1, aU2, aV2);
     //
-    //printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
-    printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
-          i, aX, aY, aZ, aU1, aV1, aU2, aV2);
+    printf("point p_%d %lf %lf %lf\n", i, aX, aY, aZ);
+    //printf("point p_%d %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf %20.15lf\n",
+       //   i, aX, aY, aZ, aU1, aV1, aU2, aV2);
   }
 }
-//modified by NIZNHY-PKV Fri Nov 25 11:58:19 2011t
+//=======================================================================
+//function : RefineVector
+//purpose  : 
+//=======================================================================
+void RefineVector(gp_Vec2d& aV2D)
+{
+  Standard_Integer k,m;
+  Standard_Real aC[2], aEps, aR1, aR2, aNum;
+  //
+  aEps=RealEpsilon();
+  aR1=1.-aEps;
+  aR2=1.+aEps;
+  //
+  aV2D.Coord(aC[0], aC[1]);
+  //
+  for (k=0; k<2; ++k) {
+    m=(k+1)%2;
+    aNum=fabs(aC[k]);
+    if (aNum>aR1 && aNum<aR2) {
+      if (aC[k]<0.) {
+       aC[k]=-1.;
+      }          
+      else {
+       aC[k]=1.;
+      }
+      aC[m]=0.;
+      break;
+    }
+  }
+  aV2D.SetCoord(aC[0], aC[1]);
+}
+//=======================================================================
+//function : FindMaxSquareDistance
+//purpose  : 
+//=======================================================================
+Standard_Real FindMaxSquareDistance (const Standard_Real aT1,
+                                    const Standard_Real aT2,
+                                    const Standard_Real aEps,
+                                    const Handle(Geom_Curve)& aC3D,
+                                    const Handle(Geom2d_Curve)& aC2D1,
+                                    const Handle(Geom2d_Curve)& aC2D2,
+                                    const Handle(GeomAdaptor_HSurface)& myHS1,
+                                    const Handle(GeomAdaptor_HSurface)& myHS2,
+                                    const TopoDS_Face& myFace1,
+                                    const TopoDS_Face& myFace2,
+                                    const Handle(BOPInt_Context)& myContext)
+{
+  Standard_Real aA, aB, aCf, aX1, aX2, aF1, aF2, aX, aF;
+  //
+  aCf=1.6180339887498948482045868343656;// =0.5*(1.+sqrt(5.));
+  aA=aT1;
+  aB=aT2;
+  aX1=aB-(aB-aA)/aCf;  
+  aF1=MaxSquareDistance(aX1, 
+                       aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+  aX2=aA+(aB-aA)/aCf;
+  aF2=MaxSquareDistance(aX2, 
+                       aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+  //
+  while(1) {
+    //
+    if (fabs(aA-aB)<aEps) {
+      aX=0.5*(aA+aB);
+      aF=MaxSquareDistance(aX, 
+                       aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+      break;
+    }
+    if (aF1<aF2){
+      aA=aX1;
+      aX1=aX2;
+      aF1=aF2;
+      aX2=aA+(aB-aA)/aCf;
+      aF2=MaxSquareDistance(aX2, 
+                           aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+      
+    }
+    else {
+      aB=aX2;
+      aX2=aX1;
+      aF2=aF1;
+      aX1=aB-(aB-aA)/aCf; 
+      aF1=MaxSquareDistance(aX1, 
+                           aC3D, aC2D1, aC2D2, myHS1, myHS2, myFace1, myFace2, myContext);
+    }
+  }
+  return aF;
+}
+//=======================================================================
+//function : MaxSquareDistance
+//purpose  : 
+//=======================================================================
+Standard_Real MaxSquareDistance (const Standard_Real aT,
+                                const Handle(Geom_Curve)& aC3D,
+                                const Handle(Geom2d_Curve)& aC2D1,
+                                const Handle(Geom2d_Curve)& aC2D2,
+                                const Handle(GeomAdaptor_HSurface) myHS1,
+                                const Handle(GeomAdaptor_HSurface) myHS2,
+                                const TopoDS_Face& aF1,
+                                const TopoDS_Face& aF2,
+                                const Handle(BOPInt_Context)& aCtx)
+{
+  Standard_Boolean bIsDone;
+  Standard_Integer i;
+  Standard_Real aU, aV, aD2Max, aD2;
+  gp_Pnt2d aP2D;
+  gp_Pnt aP, aPS;
+  //
+  aD2Max=0.;
+  //
+  aC3D->D0(aT, aP);
+  if (aC3D.IsNull()) {
+    return aD2Max;
+  }
+  //
+  for (i=0; i<2; ++i) {
+    const Handle(GeomAdaptor_HSurface)& aGHS=(!i) ? myHS1 : myHS2;
+    const TopoDS_Face &aF=(!i) ? aF1 : aF2;
+    const Handle(Geom2d_Curve)& aC2D=(!i) ? aC2D1 : aC2D2;
+    //
+    if (!aC2D.IsNull()) {
+      aC2D->D0(aT, aP2D);
+      aP2D.Coord(aU, aV);
+      aGHS->D0(aU, aV, aPS);
+      aD2=aP.SquareDistance(aPS);
+      if (aD2>aD2Max) {
+       aD2Max=aD2;
+      }
+    }
+    //
+    GeomAPI_ProjectPointOnSurf& aProjector=aCtx->ProjPS(aF);
+    //
+    aProjector.Perform(aP);
+    bIsDone=aProjector.IsDone();
+    if (bIsDone) {
+      aProjector.LowerDistanceParameters(aU, aV);
+      aGHS->D0(aU, aV, aPS);
+      aD2=aP.SquareDistance(aPS);
+      if (aD2>aD2Max) {
+       aD2Max=aD2;
+      }
+    }
+  }
+  //
+  return aD2Max;
+}
+
+//=======================================================================
+//function : CheckPCurve
+//purpose  : Checks if points of the pcurve are out of the face bounds.
+//=======================================================================
+  Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
+                               const TopoDS_Face& aFace) 
+{
+  const Standard_Integer NPoints = 23;
+  Standard_Real umin,umax,vmin,vmax;
+
+  BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
+  Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
+  Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
+  Standard_Real fp = aPC->FirstParameter();
+  Standard_Real lp = aPC->LastParameter();
+  Standard_Real step = (lp-fp)/(NPoints+1);
+
+  // adjust domain for periodic surfaces
+  TopLoc_Location aLoc;
+  Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
+  if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
+    aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
+
+  gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
+  Standard_Real u,v;
+  pnt.Coord(u,v);
+
+  if (aSurf->IsUPeriodic()) {
+    Standard_Real aPer = aSurf->UPeriod();
+    Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
+    if (u < umin+aPer*nshift) nshift--;
+    umin += aPer*nshift;
+    umax += aPer*nshift;
+  }
+  if (aSurf->IsVPeriodic()) {
+    Standard_Real aPer = aSurf->VPeriod();
+    Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
+    if (v < vmin+aPer*nshift) nshift--;
+    vmin += aPer*nshift;
+    vmax += aPer*nshift;
+  }
+
+  Standard_Integer i;
+  for (i=1; i <= NPoints; i++) {
+    Standard_Real p = fp + i * step;
+    pnt = aPC->Value(p);
+    pnt.Coord(u,v);
+    if (umin-u > tolU || u-umax > tolU ||
+        vmin-v > tolV || v-vmax > tolV)
+      return Standard_False;
+  }
+  return Standard_True;
+
+}