0024023: Revamp the OCCT Handle -- downcast (automatic)
[occt.git] / src / IntTools / IntTools_FaceFace.cxx
index ad4a922..0785460 100644 (file)
@@ -87,6 +87,7 @@
 
 #include <Geom2dAPI_InterCurveCurve.hxx>
 #include <Geom2dInt_GInter.hxx>
+#include <Geom2dAdaptor.hxx>
 #include <GeomAdaptor_Curve.hxx>
 #include <GeomAdaptor_HSurface.hxx>
 #include <GeomAdaptor_Surface.hxx>
 #include <IntTools_TopolTool.hxx>
 #include <IntTools_PntOnFace.hxx>
 #include <IntTools_PntOn2Faces.hxx>
-#include <BOPInt_Context.hxx>
+#include <IntTools_Context.hxx>
 #include <IntSurf_ListIteratorOfListOfPntOn2S.hxx>
+#include <GeomInt.hxx>
+
+#include <Approx_CurveOnSurface.hxx>
+#include <GeomAdaptor.hxx>
+#include <GeomInt_IntSS.hxx>
 
 static
   void RefineVector(gp_Vec2d& aV2D);
-#ifdef DEB_DUMPWLINE
+#ifdef OCCT_DEBUG_DUMPWLINE
 static
   void DumpWLine(const Handle(IntPatch_WLine)& aWLine);
 #endif
@@ -139,12 +145,6 @@ static
                   Standard_Real&,
                   Standard_Real&);
 
-static 
-  void BuildPCurves (Standard_Real f,Standard_Real l,Standard_Real& Tol,
-                     const Handle (Geom_Surface)& S,
-                     const Handle (Geom_Curve)&   C,
-                     Handle (Geom2d_Curve)& C2d);
-
 static 
   void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
                                 const Standard_Real theTolerance,
@@ -152,6 +152,7 @@ static
                                 Standard_Real&      theumax, 
                                 Standard_Real&      thevmin, 
                                 Standard_Real&      thevmax);
+
 static
   Standard_Boolean NotUseSurfacesForApprox
           (const TopoDS_Face& aF1,
@@ -163,13 +164,6 @@ static
 static 
   Handle(IntPatch_WLine) ComputePurgedWLine(const Handle(IntPatch_WLine)& theWLine);
 
-static 
-  Standard_Real AdjustPeriodic(const Standard_Real theParameter,
-                               const Standard_Real parmin,
-                               const Standard_Real parmax,
-                               const Standard_Real thePeriod,
-                               Standard_Real&      theOffset);
-
 static 
   Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine,
                                             const Standard_Integer                         ideb,
@@ -182,11 +176,11 @@ static
                                         const Handle(GeomAdaptor_HSurface)&            theSurface2,
                                         const TopoDS_Face&                             theFace1,
                                         const TopoDS_Face&                             theFace2,
-                                        const IntTools_LineConstructor&                theLConstructor,
+                                        const GeomInt_LineConstructor&                 theLConstructor,
                                         const Standard_Boolean                         theAvoidLConstructor,
                                         IntPatch_SequenceOfLine&                       theNewLines,
                                         Standard_Real&                                 theReachedTol3d,
-                                        const Handle(BOPInt_Context)& );
+                                        const Handle(IntTools_Context)& );
 
 static 
   Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter, 
@@ -196,7 +190,7 @@ static
                                           const Standard_Real       theOtherParameter,
                                           const Standard_Boolean    bIncreasePar,
                                           Standard_Real&            theNewParameter,
-                                          const Handle(BOPInt_Context)& );
+                                          const Handle(IntTools_Context)& );
 
 static 
   Standard_Boolean IsCurveValid(Handle(Geom2d_Curve)& thePCurve);
@@ -225,7 +219,7 @@ static
                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
                                        Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
                                        Handle(TColStd_HArray1OfReal)&       theResultRadius,
-                                       const Handle(BOPInt_Context)& );
+                                       const Handle(IntTools_Context)& );
 
 static
   Standard_Boolean FindPoint(const gp_Pnt2d&     theFirstPoint,
@@ -287,34 +281,15 @@ 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);
+  void CorrectPlaneBoundaries(Standard_Real& aUmin,
+                              Standard_Real& aUmax, 
+                              Standard_Real& aVmin, 
+                              Standard_Real& aVmax);
 
 //=======================================================================
 //function : 
@@ -336,7 +311,7 @@ IntTools_FaceFace::IntTools_FaceFace()
 //function : SetContext
 //purpose  : 
 //======================================================================= 
-void IntTools_FaceFace::SetContext(const Handle(BOPInt_Context)& aContext)
+void IntTools_FaceFace::SetContext(const Handle(IntTools_Context)& aContext)
 {
   myContext=aContext;
 }
@@ -344,7 +319,7 @@ void IntTools_FaceFace::SetContext(const Handle(BOPInt_Context)& aContext)
 //function : Context
 //purpose  : 
 //======================================================================= 
-const Handle(BOPInt_Context)& IntTools_FaceFace::Context()const
+const Handle(IntTools_Context)& IntTools_FaceFace::Context()const
 {
   return myContext;
 }
@@ -512,9 +487,6 @@ static Standard_Boolean isTreatAnalityc(const TopoDS_Face& theF1,
     return inter.IsDone();
   }
 }
-
-
-
 //=======================================================================
 //function : Perform
 //purpose  : intersect surfaces of the faces
@@ -525,7 +497,7 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   Standard_Boolean RestrictLine = Standard_False, hasCone = Standard_False;
   
   if (myContext.IsNull()) {
-    myContext=new BOPInt_Context;
+    myContext=new IntTools_Context;
   }
 
   mySeqOfCurve.Clear();
@@ -563,6 +535,10 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
         aP2S.SetValue(aU2,aV2,aU1,aV1);
       }
     }
+    //
+    Standard_Boolean anAproxTmp = myApprox1;
+    myApprox1 = myApprox2;
+    myApprox2 = anAproxTmp;
   }
 
 
@@ -583,42 +559,40 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
                                         aType2 == GeomAbs_Cone ||
                                         aType2 == GeomAbs_Torus);
 
-  if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)
-  {
+  if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane)  {
     Standard_Real umin, umax, vmin, vmax;
+    //
     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
+    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     //
     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
+    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
+    //
     Standard_Real TolAng = 1.e-8;
-
-    PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2,
-                                          mySeqOfCurve, myTangentFaces);
-
+    //
+    PerformPlanes(myHS1, myHS2, TolAng, TolTang, myApprox1, myApprox2, 
+                  mySeqOfCurve, myTangentFaces);
+    //
     myIsDone = Standard_True;
     
-    if(!myTangentFaces)
-    {
+    if(!myTangentFaces) {
       const Standard_Integer NbLinPP = mySeqOfCurve.Length();
-      if(NbLinPP)
-      {
+      if(NbLinPP) {
         Standard_Real aTolFMax;
         myTolReached3d = 1.e-7;
         aTolFMax=Max(aTolF1, aTolF2);
-        if (aTolFMax>myTolReached3d)
-        {
+        if (aTolFMax>myTolReached3d) {
           myTolReached3d=aTolFMax;
         }
-
+        //
         myTolReached2d = myTolReached3d;
 
-        if (bReverse)
-        {
+        if (bReverse) {
           Handle(Geom2d_Curve) aC2D1, aC2D2;
           const Standard_Integer aNbLin = mySeqOfCurve.Length();
-          for (Standard_Integer i = 1; i <= aNbLin; ++i)
-          {
+          for (Standard_Integer i = 1; i <= aNbLin; ++i) {
             IntTools_Curve& aIC=mySeqOfCurve(i);
             aC2D1=aIC.FirstCurve2d();
             aC2D2=aIC.SecondCurve2d();
@@ -628,24 +602,15 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
         }
       }
     }
-
     return;
   }//if(aType1==GeomAbs_Plane && aType2==GeomAbs_Plane){
 
   if ((aType1==GeomAbs_Plane) && isFace2Quad)
   {
-    Standard_Real dU, dV;
-
-    // F1
     Standard_Real umin, umax, vmin, vmax;
-    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
-
-    dU=0.1*(umax-umin);
-    dV=0.1*(vmax-vmin);
-    umin=umin-dU;
-    umax=umax+dU;
-    vmin=vmin-dV;
-    vmax=vmax+dV;
+    // F1
+    BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax); 
+    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     // F2
     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
@@ -659,21 +624,14 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   }
   else if ((aType2==GeomAbs_Plane) && isFace1Quad)
   {
-    Standard_Real dU, dV;
-
-    //F1
     Standard_Real umin, umax, vmin, vmax;
+    //F1
     BRepTools::UVBounds(myFace1, umin, umax, vmin, vmax);
     CorrectSurfaceBoundaries(myFace1, (aTolF1 + aTolF2) * 2., umin, umax, vmin, vmax);
     myHS1->ChangeSurface().Load(S1, umin, umax, vmin, vmax);
     // F2
     BRepTools::UVBounds(myFace2, umin, umax, vmin, vmax);
-    dU=0.1*(umax-umin);
-    dV=0.1*(vmax-vmin);
-    umin=umin-dU;
-    umax=umax+dU;
-    vmin=vmin-dV;
-    vmax=vmax+dV;
+    CorrectPlaneBoundaries(umin, umax, vmin, vmax);
     myHS2->ChangeSurface().Load(S2, umin, umax, vmin, vmax);
     //
     if( aType1==GeomAbs_Cone ) {
@@ -703,7 +661,7 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   {
     const Standard_Real UVMaxStep = 0.001;
     const Standard_Real Deflection = (hasCone) ? 0.085 : 0.1;
-    myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection);
+  myIntersector.SetTolerances(TolArc, TolTang, UVMaxStep, Deflection); 
   }
   
   if((myHS1->IsUClosed() && !myHS1->IsUPeriodic()) || 
@@ -715,10 +673,10 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   }
   //
   if((aType1 != GeomAbs_BSplineSurface) &&
-     (aType1 != GeomAbs_BezierSurface) &&
+      (aType1 != GeomAbs_BezierSurface)  &&
      (aType1 != GeomAbs_OtherSurface)  &&
      (aType2 != GeomAbs_BSplineSurface) &&
-     (aType2 != GeomAbs_BezierSurface) &&
+      (aType2 != GeomAbs_BezierSurface)  &&
      (aType2 != GeomAbs_OtherSurface))
   {
     RestrictLine = Standard_True;
@@ -790,6 +748,7 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     }
 
     // Points
+    Standard_Boolean bValid2D1, bValid2D2;
     Standard_Real U1,V1,U2,V2;
     IntTools_PntOnFace aPntOnF1, aPntOnF2;
     IntTools_PntOn2Faces aPntOn2Faces;
@@ -800,6 +759,19 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
       const IntSurf_PntOn2S& aISPnt=myIntersector.Point(i).PntOn2S();
       const gp_Pnt& aPnt=aISPnt.Value();
       aISPnt.Parameters(U1,V1,U2,V2);
+      //
+      // check the validity of the intersection point for the faces
+      bValid2D1 = myContext->IsPointInOnFace(myFace1, gp_Pnt2d(U1, V1));
+      if (!bValid2D1) {
+        continue;
+      }
+      //
+      bValid2D2 = myContext->IsPointInOnFace(myFace2, gp_Pnt2d(U2, V2));
+      if (!bValid2D2) {
+        continue;
+      }
+      //
+      // add the intersection point
       aPntOnF1.Init(myFace1, aPnt, U1, V1);
       aPntOnF2.Init(myFace2, aPnt, U2, V2);
       //
@@ -819,13 +791,65 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   }
 }
 
+//=======================================================================
+//function : ComputeTolerance
+//purpose  : 
+//=======================================================================
+Standard_Real IntTools_FaceFace::ComputeTolerance()
+{
+  Standard_Integer i, j, aNbLin;
+  Standard_Real aFirst, aLast, aD, aDMax, aT;
+  Handle(Geom_Surface) aS1, aS2;
+  //
+  aDMax = 0;
+  aNbLin = mySeqOfCurve.Length();
+  //
+  aS1 = myHS1->ChangeSurface().Surface();
+  aS2 = myHS2->ChangeSurface().Surface();
+  //
+  for (i = 1; i <= aNbLin; ++i)
+  {
+    const IntTools_Curve& aIC = mySeqOfCurve(i);
+    const Handle(Geom_Curve)& aC3D = aIC.Curve();
+    if (aC3D.IsNull())
+    {
+      continue;
+    }
+    //
+    aFirst = aC3D->FirstParameter();
+    aLast  = aC3D->LastParameter();
+    //
+    const Handle(Geom2d_Curve)& aC2D1 = aIC.FirstCurve2d();
+    const Handle(Geom2d_Curve)& aC2D2 = aIC.SecondCurve2d();
+    //
+    for (j = 0; j < 2; ++j)
+    {
+      const Handle(Geom2d_Curve)& aC2D = !j ? aC2D1 : aC2D2;
+      const Handle(Geom_Surface)& aS = !j ? aS1 : aS2;
+      //
+      if (!aC2D.IsNull())
+      {
+        if (IntTools_Tools::ComputeTolerance
+            (aC3D, aC2D, aS, aFirst, aLast, aD, aT))
+        {
+          if (aD > aDMax)
+          {
+            aDMax = aD;
+          }
+        }
+      }
+    }
+  }
+  //
+  return aDMax;
+}
+
 //=======================================================================
 //function :ComputeTolReached3d 
 //purpose  : 
 //=======================================================================
   void IntTools_FaceFace::ComputeTolReached3d()
 {
-  Standard_Boolean bCase1;
   Standard_Integer aNbLin, i;
   GeomAbs_SurfaceType aType1, aType2;
   //
@@ -837,9 +861,6 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   aType1=myHS1->Surface().GetType();
   aType2=myHS2->Surface().GetType();
   //
-  bCase1=((aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
-          (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion));
-  //
   if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
     if (aNbLin==2){ 
       Handle(IntPatch_Line) aIL1, aIL2;
@@ -868,45 +889,12 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     }
     //ZZ
     if (aNbLin) {// Check the distances
-      Standard_Integer  aNbP, j ;
-      Standard_Real aT1, aT2, dT, aD2, aD2Max, aEps, aT11, aT12;
-      //
-      aD2Max=0.;
-      aNbP=10;
-      aNbLin=mySeqOfCurve.Length();
+      Standard_Real aDMax;
       //
-      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);
+      aDMax = ComputeTolerance();
+      if (aDMax > 0.) {
+        myTolReached3d = aDMax;
+      }
     }// if (aNbLin) 
   }// if (aType1==GeomAbs_Cylinder && aType2==GeomAbs_Cylinder) {
   //
@@ -975,8 +963,7 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
     //
     const IntTools_Curve& aIC=mySeqOfCurve(1);
     const Handle(Geom_Curve)& aC3D=aIC.Curve();
-    const Handle(Geom_BSplineCurve)& aBS=
-      Handle(Geom_BSplineCurve)::DownCast(aC3D);
+    Handle(Geom_BSplineCurve) aBS (Handle(Geom_BSplineCurve)::DownCast(aC3D));
     if (aBS.IsNull()) {
       return;
     }
@@ -1020,151 +1007,24 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   }// if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Torus) ||
   //
   else if ((aType1==GeomAbs_SurfaceOfRevolution && aType2==GeomAbs_Cylinder) ||
-           (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder)) {
-    Standard_Integer j, aNbP;
-    Standard_Real aT, aT1, aT2, dT, aD2max, aD2;
-    //
-    aNbLin=mySeqOfCurve.Length();
-    aD2max=0.;
-    aNbP=11;
-    //
-    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()) {
-        return;
-      }
-      //
-      aT1=aBC->FirstParameter();
-      aT2=aBC->LastParameter();
-      //
-      dT=(aT2-aT1)/(aNbP-1);
-      for (j=0; j<aNbP; ++j) {
-        aT=aT1+j*dT;
-        if (j==aNbP-1) {
-          aT=aT2;
-        }
-        //
-        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) {
-    //
-    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  j, aNbP;
-    Standard_Real aT1, aT2, dT, aD2max, aD2, aEps, aT11, aT12;
-    //
-    aNbLin=mySeqOfCurve.Length();
-    aD2max=0.;
-    aNbP=10;
+           (aType2==GeomAbs_SurfaceOfRevolution && aType1==GeomAbs_Cylinder) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_Sphere) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_SurfaceOfExtrusion) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_SurfaceOfExtrusion) ||
+           (aType1==GeomAbs_Plane && aType2==GeomAbs_BSplineSurface) ||
+           (aType2==GeomAbs_Plane && aType1==GeomAbs_BSplineSurface) ||
+           !myApprox) {
     //
-    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) {
+    Standard_Real aDMax;
     //
-    aD2=myTolReached3d*myTolReached3d;
-    if (aD2max > aD2) {
-      myTolReached3d=sqrt(aD2max);
+    aDMax = ComputeTolerance();
+    if (aDMax > myTolReached3d) {
+      myTolReached3d = aDMax;
     }
-  }//else if ((aType1==GeomAbs_Plane && aType2==GeomAbs_Sphere) ...
-  else if (!myApprox || bCase1) {
-  //else if (!myApprox) {
-    Standard_Integer 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);
   }
 }
+
 //=======================================================================
 //function : MakeCurve
 //purpose  : 
@@ -1189,9 +1049,9 @@ void IntTools_FaceFace::Perform(const TopoDS_Face& aF1,
   reApprox = Standard_False;
   //
   bPCurvesOk = Standard_True;
-
-reapprox:;
-
+ reapprox:;
+  
   Tolpc = myTolApprox;
   bAvoidLineConstructor = Standard_False;
   L = myIntersector.Line(Index);
@@ -1200,8 +1060,7 @@ reapprox:;
   if(typl==IntPatch_Walking) {
     Handle(IntPatch_Line) anewL;
     //
-    const Handle(IntPatch_WLine)& aWLine=
-      Handle(IntPatch_WLine)::DownCast(L);
+    Handle(IntPatch_WLine) aWLine (Handle(IntPatch_WLine)::DownCast(L));
     //DumpWLine(aWLine);
 
     anewL = ComputePurgedWLine(aWLine);
@@ -1210,7 +1069,7 @@ reapprox:;
     }
     L = anewL;
     
-    //const Handle(IntPatch_WLine)& aWLineX = Handle(IntPatch_WLine)::DownCast(L);
+    //Handle(IntPatch_WLine) aWLineX (Handle(IntPatch_WLine)::DownCast(L));
     //DumpWLine(aWLineX);
 
     //
@@ -1229,25 +1088,36 @@ reapprox:;
       bAvoidLineConstructor = Standard_False;
     }
   }
+
+  typl=L->ArcType();
+
   //
   // Line Constructor
   if(!bAvoidLineConstructor) {
     myLConstruct.Perform(L);
     //
     bDone=myLConstruct.IsDone();
-    aNbParts=myLConstruct.NbParts();
-    if (!bDone|| !aNbParts) {
+    if(!bDone)
+    {
       return;
     }
+
+    if(typl != IntPatch_Restriction)
+    {
+      aNbParts=myLConstruct.NbParts();
+      if (aNbParts <= 0)
+      {
+        return;
+      }
+    }
   }
   // Do the Curve
-
-
-  typl=L->ArcType();
+  
+  
   switch (typl) {
-    //########################################  
-    // Line, Parabola, Hyperbola
-    //########################################  
+  //########################################  
+  // Line, Parabola, Hyperbola
+  //########################################  
   case IntPatch_Lin:
   case IntPatch_Parabola: 
   case IntPatch_Hyperbola: {
@@ -1260,7 +1130,7 @@ reapprox:;
       newc = 
         new Geom_Parabola(Handle(IntPatch_GLine)::DownCast(L)->Parabola());
     }
-
+    
     else if (typl == IntPatch_Hyperbola) {
       newc = 
         new Geom_Hyperbola (Handle(IntPatch_GLine)::DownCast(L)->Hyperbola());
@@ -1273,10 +1143,14 @@ reapprox:;
     //
     aNbParts=myLConstruct.NbParts();
     for (i=1; i<=aNbParts; i++) {
+      Standard_Boolean bFNIt, bLPIt;
+      //
       myLConstruct.Part(i, fprm, lprm);
-
-      if (!Precision::IsNegativeInfinite(fprm) && 
-          !Precision::IsPositiveInfinite(lprm)) {
+        //
+      bFNIt=Precision::IsNegativeInfinite(fprm);
+      bLPIt=Precision::IsPositiveInfinite(lprm);
+      //
+      if (!bFNIt && !bLPIt) {
         //
         IntTools_Curve aCurve;
         //
@@ -1284,7 +1158,7 @@ reapprox:;
         aCurve.SetCurve(aCT3D);
         if (typl == IntPatch_Parabola) {
           Standard_Real aTolF1, aTolF2, aTolBase;
-
+          
           aTolF1 = BRep_Tool::Tolerance(myFace1);
           aTolF2 = BRep_Tool::Tolerance(myFace2);
           aTolBase=aTolF1+aTolF2;
@@ -1294,53 +1168,54 @@ reapprox:;
         aCurve.SetCurve(new Geom_TrimmedCurve(newc, fprm, lprm));
         if(myApprox1) { 
           Handle (Geom2d_Curve) C2d;
-          BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
+          GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
+                myHS1->ChangeSurface().Surface(), newc, C2d);
           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
             myTolReached2d=Tolpc;
           }
-          //            
-          aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
-        }
-        else { 
-          Handle(Geom2d_BSplineCurve) H1;
-          //
-          aCurve.SetFirstCurve2d(H1);
-        }
-
+            //            
+            aCurve.SetFirstCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
+          }
+          else { 
+            Handle(Geom2d_BSplineCurve) H1;
+            //
+            aCurve.SetFirstCurve2d(H1);
+          }
+        //
         if(myApprox2) { 
           Handle (Geom2d_Curve) C2d;
-          BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
+          GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
+                    myHS2->ChangeSurface().Surface(), newc, C2d);
           if(Tolpc>myTolReached2d || myTolReached2d==0.) { 
             myTolReached2d=Tolpc;
           }
           //
           aCurve.SetSecondCurve2d(new Geom2d_TrimmedCurve(C2d,fprm,lprm));
-        }
+          }
         else { 
           Handle(Geom2d_BSplineCurve) H1;
           //
           aCurve.SetSecondCurve2d(H1);
         }
         mySeqOfCurve.Append(aCurve);
-      } // end of if (!Precision::IsNegativeInfinite(fprm) &&  !Precision::IsPositiveInfinite(lprm))
+      } //if (!bFNIt && !bLPIt) {
       else {
         //  on regarde si on garde
         //
-        Standard_Boolean bFNIt, bLPIt;
         Standard_Real aTestPrm, dT=100.;
-
-        bFNIt=Precision::IsNegativeInfinite(fprm);
-        bLPIt=Precision::IsPositiveInfinite(lprm);
-
+        //
         aTestPrm=0.;
-
         if (bFNIt && !bLPIt) {
           aTestPrm=lprm-dT;
         }
         else if (!bFNIt && bLPIt) {
           aTestPrm=fprm+dT;
         }
-
+        else {
+          // i.e, if (bFNIt && bLPIt)
+          aTestPrm=IntTools_Tools::IntermediatePoint(-dT, dT);
+        }
+        //
         gp_Pnt ptref(newc->Value(aTestPrm));
         //
         GeomAbs_SurfaceType typS1 = myHS1->GetType();
@@ -1350,15 +1225,14 @@ reapprox:;
             typS1 == GeomAbs_SurfaceOfRevolution ||
             typS2 == GeomAbs_SurfaceOfExtrusion ||
             typS2 == GeomAbs_OffsetSurface ||
-            typS2 == GeomAbs_SurfaceOfRevolution) 
-        {
+            typS2 == GeomAbs_SurfaceOfRevolution) {
           Handle(Geom2d_BSplineCurve) H1;
           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
           continue;
         }
 
         Standard_Real u1, v1, u2, v2, Tol;
-
+        
         Tol = Precision::Confusion();
         Parameters(myHS1, myHS2, ptref,  u1, v1, u2, v2);
         ok = (dom1->Classify(gp_Pnt2d(u1, v1), Tol) != TopAbs_OUT);
@@ -1370,13 +1244,13 @@ reapprox:;
           mySeqOfCurve.Append(IntTools_Curve(newc, H1, H1));
         }
       }
-    }// end of for (i=1; i<=myLConstruct.NbParts(); i++)
+    }// for (i=1; i<=aNbParts; i++) {
   }// case IntPatch_Lin:  case IntPatch_Parabola:  case IntPatch_Hyperbola:
-  break;
+    break;
 
-    //########################################  
-    // Circle and Ellipse
-    //########################################  
+  //########################################  
+  // Circle and Ellipse
+  //########################################  
   case IntPatch_Circle: 
   case IntPatch_Ellipse: {
 
@@ -1396,7 +1270,7 @@ reapprox:;
     //
     Standard_Real aPeriod, aNul;
     TColStd_SequenceOfReal aSeqFprm,  aSeqLprm;
-
+    
     aNul=0.;
     aPeriod=M_PI+M_PI;
 
@@ -1459,7 +1333,7 @@ reapprox:;
         aSeqLprm.Append(lprm);
       }
     }
-
+    
     //
     aNbParts=aSeqFprm.Length();
     for (i=1; i<=aNbParts; i++) {
@@ -1479,7 +1353,8 @@ reapprox:;
         if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {//// 
           if(myApprox1) { 
             Handle (Geom2d_Curve) C2d;
-            BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
+            GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
+                        myHS1->ChangeSurface().Surface(), newc, C2d);
             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
               myTolReached2d=Tolpc;
             }
@@ -1494,7 +1369,8 @@ reapprox:;
 
           if(myApprox2) { 
             Handle (Geom2d_Curve) C2d;
-            BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
+            GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
+                    myHS2->ChangeSurface().Surface(),newc,C2d);
             if(Tolpc>myTolReached2d || myTolReached2d==0) { 
               myTolReached2d=Tolpc;
             }
@@ -1506,31 +1382,32 @@ reapprox:;
             aCurve.SetSecondCurve2d(H1);
           }
         }
-
+        
         else { 
           Handle(Geom2d_BSplineCurve) H1;
           aCurve.SetFirstCurve2d(H1);
           aCurve.SetSecondCurve2d(H1);
         }
         mySeqOfCurve.Append(aCurve);
-        //==============================================        
+          //==============================================        
       } //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.*M_PI) < RealEpsilon()) {
+//           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);
             fprm=aTC3D->FirstParameter();
             lprm=aTC3D->LastParameter ();
-
+            
             if(myApprox1) { 
               Handle (Geom2d_Curve) C2d;
-              BuildPCurves(fprm,lprm,Tolpc,myHS1->ChangeSurface().Surface(),newc,C2d);
+              GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
+                    myHS1->ChangeSurface().Surface(),newc,C2d);
               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
                 myTolReached2d=Tolpc;
               }
@@ -1544,7 +1421,8 @@ reapprox:;
 
             if(myApprox2) { 
               Handle (Geom2d_Curve) C2d;
-              BuildPCurves(fprm,lprm,Tolpc,myHS2->ChangeSurface().Surface(),newc,C2d);
+              GeomInt_IntSS::BuildPCurves(fprm,lprm,Tolpc,
+                    myHS2->ChangeSurface().Surface(),newc,C2d);
               if(Tolpc>myTolReached2d || myTolReached2d==0) { 
                 myTolReached2d=Tolpc;
               }
@@ -1578,10 +1456,11 @@ reapprox:;
             aCurve.SetCurve(newc);
             //==============================================
             if (typl == IntPatch_Circle || typl == IntPatch_Ellipse) {
-
+              
               if(myApprox1) { 
                 Handle (Geom2d_Curve) C2d;
-                BuildPCurves(fprm, lprm, Tolpc, myHS1->ChangeSurface().Surface(), newc, C2d);
+                GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc, 
+                        myHS1->ChangeSurface().Surface(), newc, C2d);
                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
                   myTolReached2d=Tolpc;
                 }
@@ -1592,23 +1471,24 @@ reapprox:;
                 Handle(Geom2d_BSplineCurve) H1;
                 aCurve.SetFirstCurve2d(H1);
               }
-
+                
               if(myApprox2) { 
                 Handle (Geom2d_Curve) C2d;
-                BuildPCurves(fprm, lprm, Tolpc,myHS2->ChangeSurface().Surface(), newc, C2d);
+                GeomInt_IntSS::BuildPCurves(fprm, lprm, Tolpc,
+                        myHS2->ChangeSurface().Surface(), newc, C2d);
                 if(Tolpc>myTolReached2d || myTolReached2d==0) { 
                   myTolReached2d=Tolpc;
                 }
                 //                 
                 aCurve.SetSecondCurve2d(C2d);
               }
-
+                
               else { 
                 Handle(Geom2d_BSplineCurve) H1;
                 aCurve.SetSecondCurve2d(H1);
               }
             }//  end of if (typl == IntPatch_Circle || typl == IntPatch_Ellipse)
-
+             
             else { 
               Handle(Geom2d_BSplineCurve) H1;
               //         
@@ -1620,62 +1500,62 @@ reapprox:;
             mySeqOfCurve.Append(aCurve);
             break;
 
-          }//  end of if (ok) {
-        }//  end of for (Standard_Integer j=0; j<=17; j++)
-      }//  end of else { on regarde si on garde
-    }// for (i=1; i<=myLConstruct.NbParts(); i++)
-  }// IntPatch_Circle: IntPatch_Ellipse:
-  break;
-
+            }//  end of if (ok) {
+          }//  end of for (Standard_Integer j=0; j<=17; j++)
+        }//  end of else { on regarde si on garde
+      }// for (i=1; i<=myLConstruct.NbParts(); i++)
+    }// IntPatch_Circle: IntPatch_Ellipse:
+    break;
+    
   case IntPatch_Analytic: {
     IntSurf_Quadric quad1,quad2;
     GeomAbs_SurfaceType typs = myHS1->Surface().GetType();
-
+    
     switch (typs) {
-    case GeomAbs_Plane:
-      quad1.SetValue(myHS1->Surface().Plane());
-      break;
-    case GeomAbs_Cylinder:
-      quad1.SetValue(myHS1->Surface().Cylinder());
-      break;
-    case GeomAbs_Cone:
-      quad1.SetValue(myHS1->Surface().Cone());
-      break;
-    case GeomAbs_Sphere:
-      quad1.SetValue(myHS1->Surface().Sphere());
-      break;
+      case GeomAbs_Plane:
+        quad1.SetValue(myHS1->Surface().Plane());
+        break;
+      case GeomAbs_Cylinder:
+        quad1.SetValue(myHS1->Surface().Cylinder());
+        break;
+      case GeomAbs_Cone:
+        quad1.SetValue(myHS1->Surface().Cone());
+        break;
+      case GeomAbs_Sphere:
+        quad1.SetValue(myHS1->Surface().Sphere());
+        break;
     case GeomAbs_Torus:
       quad1.SetValue(myHS1->Surface().Torus());
       break;
-    default:
-      Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
-    }
-
+      default:
+        Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 1");
+      }
+      
     typs = myHS2->Surface().GetType();
-
+    
     switch (typs) {
-    case GeomAbs_Plane:
-      quad2.SetValue(myHS2->Surface().Plane());
-      break;
-    case GeomAbs_Cylinder:
-      quad2.SetValue(myHS2->Surface().Cylinder());
-      break;
-    case GeomAbs_Cone:
-      quad2.SetValue(myHS2->Surface().Cone());
-      break;
-    case GeomAbs_Sphere:
-      quad2.SetValue(myHS2->Surface().Sphere());
-      break;
+      case GeomAbs_Plane:
+        quad2.SetValue(myHS2->Surface().Plane());
+        break;
+      case GeomAbs_Cylinder:
+        quad2.SetValue(myHS2->Surface().Cylinder());
+        break;
+      case GeomAbs_Cone:
+        quad2.SetValue(myHS2->Surface().Cone());
+        break;
+      case GeomAbs_Sphere:
+        quad2.SetValue(myHS2->Surface().Sphere());
+        break;
     case GeomAbs_Torus:
       quad2.SetValue(myHS2->Surface().Torus());
       break;
-    default:
-      Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
-    }
+      default:
+        Standard_ConstructionError::Raise("GeomInt_IntSS::MakeCurve 2");
+      }
     //
     //=========
     IntPatch_ALineToWLine convert (quad1, quad2);
-
+      
     if (!myApprox) {
       aNbParts=myLConstruct.NbParts();
       for (i=1; i<=aNbParts; i++) {
@@ -1689,7 +1569,7 @@ reapprox:;
         if(myApprox1) {
           H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
         }
-
+        
         if(myApprox2) {
           H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
         }
@@ -1704,7 +1584,7 @@ reapprox:;
       Standard_Real tol2d = myTolApprox;
       //         
       theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True);
-
+      
       aNbParts=myLConstruct.NbParts();
       for (i=1; i<=aNbParts; i++) {
         myLConstruct.Part(i, fprm, lprm);
@@ -1721,7 +1601,7 @@ reapprox:;
           if(myApprox1) {
             H1 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_True);
           }
-
+          
           if(myApprox2) {
             H2 = MakeBSpline2d(WL, 1, WL->NbPnts(), Standard_False);
           }
@@ -1735,7 +1615,7 @@ reapprox:;
               myTolReached2d = theapp3d.TolReached2d();
             }
           }
-
+          
           if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0) { 
             myTolReached3d = theapp3d.TolReached3d();
           }
@@ -1745,27 +1625,27 @@ reapprox:;
           for (j=1; j<=aNbMultiCurves; j++) {
             const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
             nbpoles = mbspc.NbPoles();
-
+            
             TColgp_Array1OfPnt tpoles(1, nbpoles);
             mbspc.Curve(1, tpoles);
             Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
-              mbspc.Knots(),
-              mbspc.Multiplicities(),
-              mbspc.Degree());
-
+                                                               mbspc.Knots(),
+                                                               mbspc.Multiplicities(),
+                                                               mbspc.Degree());
+            
             GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
             Check.FixTangent(Standard_True,Standard_True);
             // 
             IntTools_Curve aCurve;
             aCurve.SetCurve(BS);
-
+            
             if(myApprox1) { 
               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
               mbspc.Curve(2,tpoles2d);
               Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
-                mbspc.Knots(),
-                mbspc.Multiplicities(),
-                mbspc.Degree());
+                                                                      mbspc.Knots(),
+                                                                      mbspc.Multiplicities(),
+                                                                      mbspc.Degree());
 
               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
               newCheck.FixTangent(Standard_True,Standard_True);
@@ -1776,17 +1656,17 @@ reapprox:;
               Handle(Geom2d_BSplineCurve) H1;
               aCurve.SetFirstCurve2d(H1);
             }
-
+            
             if(myApprox2) { 
               TColgp_Array1OfPnt2d tpoles2d(1, nbpoles);
               Standard_Integer TwoOrThree;
               TwoOrThree=myApprox1 ? 3 : 2;
               mbspc.Curve(TwoOrThree, tpoles2d);
               Handle(Geom2d_BSplineCurve) BS2 =new Geom2d_BSplineCurve(tpoles2d,
-                mbspc.Knots(),
-                mbspc.Multiplicities(),
-                mbspc.Degree());
-
+                                                                       mbspc.Knots(),
+                                                                       mbspc.Multiplicities(),
+                                                                       mbspc.Degree());
+                
               GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
               newCheck.FixTangent(Standard_True,Standard_True);
               //         
@@ -1804,9 +1684,9 @@ reapprox:;
       }// for (i=1; i<=aNbParts; i++) {
     }// else { // myApprox=TRUE
   }// case IntPatch_Analytic:
-  break;
+    break;
 
-  case IntPatch_Walking: {
+  case IntPatch_Walking:{
     Handle(IntPatch_WLine) WL = 
       Handle(IntPatch_WLine)::DownCast(L);
     //
@@ -1880,7 +1760,7 @@ reapprox:;
         //         
         tol2d = myTolApprox;
       }
-
+        
       if(myHS1 == myHS2) { 
         theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_False, aParType);
         rejectSurface = Standard_True;
@@ -1955,10 +1835,10 @@ reapprox:;
           //
           if (myHS1 != myHS2){
             if ((typs1==GeomAbs_BezierSurface || typs1==GeomAbs_BSplineSurface) &&
-              (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
-
+                (typs2==GeomAbs_BezierSurface || typs2==GeomAbs_BSplineSurface)) {
+             
               theapp3d.SetParameters(myTolApprox, tol2d, 4, 8, 0, Standard_True, aParType);
-
+              
               Standard_Boolean bUseSurfaces;
               bUseSurfaces=NotUseSurfacesForApprox(myFace1, myFace2, WL, ifprm,  ilprm);
               if (bUseSurfaces) {
@@ -1972,7 +1852,7 @@ reapprox:;
           //
           theapp3d.Perform(myHS1,myHS2,WL,Standard_True,anApprox1,anApprox2,ifprm,ilprm);
         }
-        //
+          //           
         if (!theapp3d.IsDone()) {
           Handle(Geom2d_BSplineCurve) H1;
           Handle(Geom2d_BSplineCurve) H2;
@@ -1990,7 +1870,7 @@ reapprox:;
           IntTools_Curve aIC(aBSp, H1, H2);
           mySeqOfCurve.Append(aIC);
         }
-
+        
         else {
           if(myApprox1 || myApprox2 || (typs1==GeomAbs_Plane || typs2==GeomAbs_Plane)) { 
             if( theapp3d.TolReached2d()>myTolReached2d || myTolReached2d==0.) { 
@@ -2010,33 +1890,33 @@ reapprox:;
           else  if( theapp3d.TolReached3d()>myTolReached3d || myTolReached3d==0.) { 
             myTolReached3d = theapp3d.TolReached3d();
           }
-
+          
           Standard_Integer aNbMultiCurves, nbpoles;
           aNbMultiCurves=theapp3d.NbMultiCurves(); 
           for (j=1; j<=aNbMultiCurves; j++) {
             if(typs1 == GeomAbs_Plane) {
               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
               nbpoles = mbspc.NbPoles();
-
+              
               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
               TColgp_Array1OfPnt   tpoles(1,nbpoles);
-
+              
               mbspc.Curve(1,tpoles2d);
               const gp_Pln&  Pln = myHS1->Surface().Plane();
               //
               Standard_Integer ik; 
               for(ik = 1; ik<= nbpoles; ik++) { 
                 tpoles.SetValue(ik,
-                  ElSLib::Value(tpoles2d.Value(ik).X(),
-                  tpoles2d.Value(ik).Y(),
-                  Pln));
+                                ElSLib::Value(tpoles2d.Value(ik).X(),
+                                              tpoles2d.Value(ik).Y(),
+                                              Pln));
               }
               //
               Handle(Geom_BSplineCurve) BS = 
                 new Geom_BSplineCurve(tpoles,
-                mbspc.Knots(),
-                mbspc.Multiplicities(),
-                mbspc.Degree());
+                                      mbspc.Knots(),
+                                      mbspc.Multiplicities(),
+                                      mbspc.Degree());
               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
               Check.FixTangent(Standard_True, Standard_True);
               //         
@@ -2046,9 +1926,9 @@ reapprox:;
               if(myApprox1) { 
                 Handle(Geom2d_BSplineCurve) BS1 = 
                   new Geom2d_BSplineCurve(tpoles2d,
-                  mbspc.Knots(),
-                  mbspc.Multiplicities(),
-                  mbspc.Degree());
+                                          mbspc.Knots(),
+                                          mbspc.Multiplicities(),
+                                          mbspc.Degree());
                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
                 Check1.FixTangent(Standard_True,Standard_True);
                 //
@@ -2070,14 +1950,14 @@ reapprox:;
 
               if(myApprox2) { 
                 mbspc.Curve(2, tpoles2d);
-
+                
                 Handle(Geom2d_BSplineCurve) BS2 = new Geom2d_BSplineCurve(tpoles2d,
-                  mbspc.Knots(),
-                  mbspc.Multiplicities(),
-                  mbspc.Degree());
+                                                                          mbspc.Knots(),
+                                                                          mbspc.Multiplicities(),
+                                                                          mbspc.Degree());
                 GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
                 newCheck.FixTangent(Standard_True,Standard_True);
-
+                
                 // ###########################################
                 if(!rejectSurface && !reApprox) {
                   Standard_Boolean isValid = IsCurveValid(BS2);
@@ -2093,16 +1973,16 @@ reapprox:;
               else { 
                 Handle(Geom2d_BSplineCurve) H2;
                 //                 
-                aCurve.SetSecondCurve2d(H2);
+                  aCurve.SetSecondCurve2d(H2);
               }
               //
               mySeqOfCurve.Append(aCurve);
             }//if(typs1 == GeomAbs_Plane) {
-
+            
             else if(typs2 == GeomAbs_Plane) { 
               const AppParCurves_MultiBSpCurve& mbspc = theapp3d.Value(j);
               nbpoles = mbspc.NbPoles();
-
+              
               TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
               TColgp_Array1OfPnt   tpoles(1,nbpoles);
               mbspc.Curve((myApprox1==Standard_True)? 2 : 1,tpoles2d);
@@ -2111,16 +1991,16 @@ reapprox:;
               Standard_Integer ik; 
               for(ik = 1; ik<= nbpoles; ik++) { 
                 tpoles.SetValue(ik,
-                  ElSLib::Value(tpoles2d.Value(ik).X(),
-                  tpoles2d.Value(ik).Y(),
-                  Pln));
-
+                                ElSLib::Value(tpoles2d.Value(ik).X(),
+                                              tpoles2d.Value(ik).Y(),
+                                              Pln));
+                
               }
               //
               Handle(Geom_BSplineCurve) BS=new Geom_BSplineCurve(tpoles,
-                mbspc.Knots(),
-                mbspc.Multiplicities(),
-                mbspc.Degree());
+                                                                 mbspc.Knots(),
+                                                                 mbspc.Multiplicities(),
+                                                                 mbspc.Degree());
               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
               Check.FixTangent(Standard_True,Standard_True);
               //         
@@ -2129,9 +2009,9 @@ reapprox:;
 
               if(myApprox2) {
                 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(tpoles2d,
-                  mbspc.Knots(),
-                  mbspc.Multiplicities(),
-                  mbspc.Degree());
+                                                                        mbspc.Knots(),
+                                                                        mbspc.Multiplicities(),
+                                                                        mbspc.Degree());
                 GeomLib_Check2dBSplineCurve Check1(BS1,TOLCHECK,TOLANGCHECK);
                 Check1.FixTangent(Standard_True,Standard_True);
                 //         
@@ -2151,13 +2031,13 @@ reapprox:;
                 Handle(Geom2d_BSplineCurve) H2;
                 aCurve.SetSecondCurve2d(H2);
               }
-
+              
               if(myApprox1) { 
                 mbspc.Curve(1,tpoles2d);
                 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(tpoles2d,
-                  mbspc.Knots(),
-                  mbspc.Multiplicities(),
-                  mbspc.Degree());
+                                                                        mbspc.Knots(),
+                                                                        mbspc.Multiplicities(),
+                                                                        mbspc.Degree());
                 GeomLib_Check2dBSplineCurve Check2(BS2,TOLCHECK,TOLANGCHECK);
                 Check2.FixTangent(Standard_True,Standard_True);
                 //
@@ -2186,12 +2066,12 @@ reapprox:;
                 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);
@@ -2224,12 +2104,12 @@ reapprox:;
               TColgp_Array1OfPnt tpoles(1,nbpoles);
               mbspc.Curve(1,tpoles);
               BS=new Geom_BSplineCurve(tpoles,
-                                       mbspc.Knots(),
-                                       mbspc.Multiplicities(),
-                                       mbspc.Degree());
+                                                                 mbspc.Knots(),
+                                                                 mbspc.Multiplicities(),
+                                                                 mbspc.Degree());
               GeomLib_CheckBSplineCurve Check(BS,TOLCHECK,TOLANGCHECK);
               Check.FixTangent(Standard_True,Standard_True);
-              //
+              //                 
               aCurve.SetCurve(BS);
               aCurve.SetFirstCurve2d(aH2D);
               aCurve.SetSecondCurve2d(aH2D);
@@ -2241,12 +2121,12 @@ reapprox:;
                   mbspc.Curve(2,tpoles2d);
                   //
                   BS1=new Geom2d_BSplineCurve(tpoles2d,
-                                              mbspc.Knots(),
-                                              mbspc.Multiplicities(),
-                                              mbspc.Degree());
+                                                                        mbspc.Knots(),
+                                                                        mbspc.Multiplicities(),
+                                                                        mbspc.Degree());
                   GeomLib_Check2dBSplineCurve newCheck(BS1,TOLCHECK,TOLANGCHECK);
                   newCheck.FixTangent(Standard_True,Standard_True);
-                  // 
+                  //         
                   if (!reApprox) {
                     bIsValid1=CheckPCurve(BS1, myFace1);
                   }
@@ -2260,24 +2140,25 @@ reapprox:;
 
                   Handle(Geom2d_Curve) C2d;
                   Standard_Real aTol = myTolApprox;
-                  BuildPCurves(fprm, lprm, aTol, myHS1->ChangeSurface().Surface(), BS, C2d);
+                  GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
+                            myHS1->ChangeSurface().Surface(), BS, C2d);
                   BS1 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
                   aCurve.SetFirstCurve2d(BS1);
                 }
               } // if(myApprox1) { 
-              //
+                //                 
               if(myApprox2) { 
                 if(anApprox2) {
                   Handle(Geom2d_BSplineCurve) BS2;
                   TColgp_Array1OfPnt2d tpoles2d(1,nbpoles);
                   mbspc.Curve((myApprox1==Standard_True)? 3 : 2,tpoles2d);
                   BS2=new Geom2d_BSplineCurve(tpoles2d,
-                                              mbspc.Knots(),
-                                              mbspc.Multiplicities(),
-                                              mbspc.Degree());
+                                                                        mbspc.Knots(),
+                                                                        mbspc.Multiplicities(),
+                                                                        mbspc.Degree());
                   GeomLib_Check2dBSplineCurve newCheck(BS2,TOLCHECK,TOLANGCHECK);
                   newCheck.FixTangent(Standard_True,Standard_True);
-                  // 
+                //                 
                   if (!reApprox) {
                     bIsValid2=CheckPCurve(BS2, myFace2);        
                   }
@@ -2290,7 +2171,8 @@ reapprox:;
 
                   Handle(Geom2d_Curve) C2d;
                   Standard_Real aTol = myTolApprox;
-                  BuildPCurves(fprm, lprm, aTol, myHS2->ChangeSurface().Surface(), BS, C2d);
+                  GeomInt_IntSS::BuildPCurves(fprm, lprm, aTol,
+                            myHS2->ChangeSurface().Surface(), BS, C2d);
                   BS2 = Handle(Geom2d_BSplineCurve)::DownCast(C2d);
                   aCurve.SetSecondCurve2d(BS2);
                 }
@@ -2301,7 +2183,7 @@ reapprox:;
                 reApprox = Standard_True;
                 goto reapprox;
               }
-              //
+                //                 
               mySeqOfCurve.Append(aCurve);
             }
           }
@@ -2309,118 +2191,117 @@ reapprox:;
       }
     }// else { // X
   }// case IntPatch_Walking:{
-  break;
-
-  case IntPatch_Restriction: 
-    break;
-  default:
     break;
+    
+  case IntPatch_Restriction: 
+    {
+      GeomAbs_SurfaceType typS1 = myHS1->Surface().GetType();
+      GeomAbs_SurfaceType typS2 = myHS2->Surface().GetType();
+      Standard_Boolean isAnalS1 = Standard_False;
+      switch (typS1)
+      {
+      case GeomAbs_Plane:
+      case GeomAbs_Cylinder:
+      case GeomAbs_Sphere:
+      case GeomAbs_Cone: 
+      case GeomAbs_Torus: isAnalS1 = Standard_True; break;
+      default: break;
+      }
 
-  }
-}
+      Standard_Integer isAnalS2 = Standard_False;
+      switch (typS2)
+      {
+      case GeomAbs_Plane:
+      case GeomAbs_Cylinder:
+      case GeomAbs_Sphere:
+      case GeomAbs_Cone: 
+      case GeomAbs_Torus: isAnalS2 = Standard_True; break;
+      default: break;
+      }
 
-//=======================================================================
-//function : BuildPCurves
-//purpose  : 
-//=======================================================================
- void BuildPCurves (Standard_Real f,
-                    Standard_Real l,
-                    Standard_Real& Tol,
-                    const Handle (Geom_Surface)& S,
-                    const Handle (Geom_Curve)&   C,
-                    Handle (Geom2d_Curve)& C2d)
-{
+      Handle(IntPatch_RLine) RL = 
+        Handle(IntPatch_RLine)::DownCast(L);
+      Handle(Geom_Curve) aC3d;
+      Handle(Geom2d_Curve) aC2d1, aC2d2;
+      Standard_Real aTolReached;
+      GeomInt_IntSS::TreatRLine(RL, myHS1, myHS2, aC3d,
+                                  aC2d1, aC2d2, aTolReached);
 
-  Standard_Real umin,umax,vmin,vmax;
-  // 
+      if(aC3d.IsNull())
+        break;
 
-  if (C2d.IsNull()) {
+      Bnd_Box2d aBox1, aBox2;
 
-    // in class ProjLib_Function the range of parameters is shrank by 1.e-09
-    if((l - f) > 2.e-09) {
-      C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
-      //
-      if (C2d.IsNull()) {
-        // proj. a circle that goes through the pole on a sphere to the sphere     
-        Tol=Tol+1.e-7;
-        C2d = GeomProjLib::Curve2d(C,f,l,S,Tol);
-      }
-    }
-    else {
-      if((l - f) > Epsilon(Abs(f))) {
-        GeomAPI_ProjectPointOnSurf aProjector1, aProjector2;
-        gp_Pnt P1 = C->Value(f);
-        gp_Pnt P2 = C->Value(l);
-        aProjector1.Init(P1, S);
-        aProjector2.Init(P2, S);
-
-        if(aProjector1.IsDone() && aProjector2.IsDone()) {
-          Standard_Real U=0., V=0.;
-          aProjector1.LowerDistanceParameters(U, V);
-          gp_Pnt2d p1(U, V);
-
-          aProjector2.LowerDistanceParameters(U, V);
-          gp_Pnt2d p2(U, V);
-
-          if(p1.Distance(p2) > gp::Resolution()) {
-            TColgp_Array1OfPnt2d poles(1,2);
-            TColStd_Array1OfReal knots(1,2);
-            TColStd_Array1OfInteger mults(1,2);
-            poles(1) = p1;
-            poles(2) = p2;
-            knots(1) = f;
-            knots(2) = l;
-            mults(1) = mults(2) = 2;
-
-            C2d = new Geom2d_BSplineCurve(poles,knots,mults,1);
-
-            // compute reached tolerance.begin
-            gp_Pnt PMid = C->Value((f + l) * 0.5);
-            aProjector1.Perform(PMid);
-
-            if(aProjector1.IsDone()) {
-              aProjector1.LowerDistanceParameters(U, V);
-              gp_Pnt2d pmidproj(U, V);
-              gp_Pnt2d pmidcurve2d = C2d->Value((f + l) * 0.5);
-              Standard_Real adist = pmidcurve2d.Distance(pmidproj);
-              Tol = (adist > Tol) ? adist : Tol;
-            }
-            // compute reached tolerance.end
-          }
+      const Standard_Real aU1f = myHS1->FirstUParameter(),
+                          aV1f = myHS1->FirstVParameter(),
+                          aU1l = myHS1->LastUParameter(),
+                          aV1l = myHS1->LastVParameter();
+      const Standard_Real aU2f = myHS2->FirstUParameter(),
+                          aV2f = myHS2->FirstVParameter(),
+                          aU2l = myHS2->LastUParameter(),
+                          aV2l = myHS2->LastVParameter();
+
+      aBox1.Add(gp_Pnt2d(aU1f, aV1f));
+      aBox1.Add(gp_Pnt2d(aU1l, aV1l));
+      aBox2.Add(gp_Pnt2d(aU2f, aV2f));
+      aBox2.Add(gp_Pnt2d(aU2l, aV2l));
+
+      GeomInt_VectorOfReal anArrayOfParameters;
+        
+      //We consider here that the intersection line is same-parameter-line
+      anArrayOfParameters.Append(aC3d->FirstParameter());
+      anArrayOfParameters.Append(aC3d->LastParameter());
+
+      GeomInt_IntSS::
+        TrimILineOnSurfBoundaries(aC2d1, aC2d2, aBox1, aBox2, anArrayOfParameters);
+
+      const Standard_Integer aNbIntersSolutionsm1 = anArrayOfParameters.Length() - 1;
+
+      //Trim RLine found.
+      for(Standard_Integer anInd = 0; anInd < aNbIntersSolutionsm1; anInd++)
+      {
+        const Standard_Real aParF = anArrayOfParameters(anInd),
+                            aParL = anArrayOfParameters(anInd+1);
+
+        if((aParL - aParF) <= Precision::PConfusion())
+          continue;
+
+        const Standard_Real aPar = 0.5*(aParF + aParL);
+        gp_Pnt2d aPt;
+
+        Handle(Geom2d_Curve) aCurv2d1, aCurv2d2;
+        if(!aC2d1.IsNull())
+        {
+          aC2d1->D0(aPar, aPt);
+
+          if(aBox1.IsOut(aPt))
+            continue;
+
+          if(myApprox1)
+            aCurv2d1 = new Geom2d_TrimmedCurve(aC2d1, aParF, aParL);
         }
-      }
-    }
-    //
-    S->Bounds(umin, umax, vmin, vmax);
 
-    if (S->IsUPeriodic() && !C2d.IsNull()) {
-      // Recadre dans le domaine UV de la face
-      Standard_Real period, U0, du, aEps; 
-      
-      du =0.0;
-      aEps=Precision::PConfusion();
-      period = S->UPeriod();
-      gp_Pnt2d Pf = C2d->Value(f);
-      U0=Pf.X();
-      //
-      gp_Pnt2d Pl = C2d->Value(l);
-      
-      U0 = Min(Pl.X(), U0);
-//       while(U0-umin<aEps) { 
-      while(U0-umin<-aEps) { 
-        U0+=period;
-        du+=period;
-      }
-      //
-      while(U0-umax>aEps) { 
-        U0-=period;
-        du-=period;
-      }
-      if (du != 0) {
-        gp_Vec2d T1(du,0.);
-        C2d->Translate(T1);
+        if(!aC2d2.IsNull())
+        {
+          aC2d2->D0(aPar, aPt);
+
+          if(aBox2.IsOut(aPt))
+            continue;
+
+          if(myApprox2)
+            aCurv2d2 = new Geom2d_TrimmedCurve(aC2d2, aParF, aParL);
+        }
+
+        Handle(Geom_Curve) aCurv3d = new Geom_TrimmedCurve(aC3d, aParF, aParL);
+
+        IntTools_Curve aIC(aCurv3d, aCurv2d1, aCurv2d2);
+        mySeqOfCurve.Append(aIC);
       }
     }
+    break;
+  default:
+    break;
+
   }
 }
 
@@ -2616,7 +2497,7 @@ Handle(Geom2d_BSplineCurve) MakeBSpline2d(const Handle(IntPatch_WLine)& theWLine
 //function : CorrectSurfaceBoundaries
 //purpose  : 
 //=======================================================================
-void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
+ void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
                               const Standard_Real theTolerance,
                               Standard_Real&      theumin,
                               Standard_Real&      theumax, 
@@ -2669,26 +2550,30 @@ void CorrectSurfaceBoundaries(const TopoDS_Face&  theFace,
   //
   if(!isuperiodic && enlarge) {
 
-    if((theumin - uinf) > delta )
+    if(!Precision::IsInfinite(theumin) &&
+        ((theumin - uinf) > delta))
       theumin -= delta;
     else {
       theumin = uinf;
     }
 
-    if((usup - theumax) > delta )
+    if(!Precision::IsInfinite(theumax) &&
+        ((usup - theumax) > delta))
       theumax += delta;
     else
       theumax = usup;
   }
   //
   if(!isvperiodic && enlarge) {
-    if((thevmin - vinf) > delta ) {
+    if(!Precision::IsInfinite(thevmin) &&
+        ((thevmin - vinf) > delta)) {
       thevmin -= delta;
     }
     else { 
       thevmin = vinf;
     }
-    if((vsup - thevmax) > delta ) {
+    if(!Precision::IsInfinite(thevmax) &&
+        ((vsup - thevmax) > delta)) {
       thevmax += delta;
     }
     else {
@@ -3118,31 +3003,6 @@ void TolR3d(const TopoDS_Face& aF1,
   }
 }
 //=======================================================================
-//function : AdjustPeriodic
-//purpose  : 
-//=======================================================================
-Standard_Real AdjustPeriodic(const Standard_Real theParameter,
-                             const Standard_Real parmin,
-                             const Standard_Real parmax,
-                             const Standard_Real thePeriod,
-                             Standard_Real&      theOffset) 
-{
-  Standard_Real aresult;
-  //
-  theOffset = 0.;
-  aresult = theParameter;
-  while(aresult < parmin) {
-    aresult += thePeriod;
-    theOffset += thePeriod;
-  }
-
-  while(aresult > parmax) {
-    aresult -= thePeriod;
-    theOffset -= thePeriod;
-  }
-  return aresult;
-}
-//=======================================================================
 //function : IsPointOnBoundary
 //purpose  : 
 //=======================================================================
@@ -3387,7 +3247,7 @@ Standard_Integer ComputeTangentZones( const Handle(GeomAdaptor_HSurface)& theSur
                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS1,
                                      Handle(TColgp_HArray1OfPnt2d)&       theResultOnS2,
                                      Handle(TColStd_HArray1OfReal)&       theResultRadius,
-                                     const Handle(BOPInt_Context)& aContext)
+                                     const Handle(IntTools_Context)& aContext)
 {
   Standard_Integer aResult = 0;
   if ( !CheckTangentZonesExist( theSurface1, theSurface2 ) )
@@ -3595,11 +3455,11 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
                                       const Handle(GeomAdaptor_HSurface)&            theSurface2,
                                       const TopoDS_Face&                             theFace1,
                                       const TopoDS_Face&                             theFace2,
-                                      const IntTools_LineConstructor&                theLConstructor,
+                                      const GeomInt_LineConstructor&                 theLConstructor,
                                       const Standard_Boolean                         theAvoidLConstructor,
                                       IntPatch_SequenceOfLine&                       theNewLines,
                                       Standard_Real&                                 theReachedTol3d,
-                                      const Handle(BOPInt_Context)& aContext) 
+                                      const Handle(IntTools_Context)& aContext) 
 {
 
   Standard_Boolean bRet, bAvoidLineConstructor;
@@ -3683,12 +3543,12 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
           aParameter=V;
         }
         
-        anoffset = 0.;
-        anAdjustPar = AdjustPeriodic(aParameter
-                                     alowerboundary, 
-                                     aupperboundary, 
-                                     aPeriod, 
-                                     anoffset);
+        GeomInt::AdjustPeriodic(aParameter, 
+                                       alowerboundary
+                                       aupperboundary, 
+                                       aPeriod,
+                                       anAdjustPar,
+                                       anoffset);
         //
         bIsOnFirstBoundary = Standard_True;// ?
         bIsPointOnBoundary=
@@ -3760,9 +3620,6 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
       continue;
     }
     const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
-    if(aListOfIndex.Extent() < 2) {
-      continue;
-    }
     TColStd_ListOfInteger aListOfFLIndex;
 
     for(j = 0; j < 2; j++) {
@@ -3778,7 +3635,15 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
       const IntSurf_PntOn2S& aPoint = theWLine->Point(anIndex);
 
       IntSurf_PntOn2S aNewP = aPoint;
-      
+      if(aListOfIndex.Extent() < 2) {
+        aSeqOfPntOn2S->Add(aNewP);
+        aListOfFLIndex.Append(aSeqOfPntOn2S->NbPoints());
+        continue;
+      }
+      //
+      Standard_Integer iFirst = aListOfIndex.First();
+      Standard_Integer iLast  = aListOfIndex.Last();
+      //
       for(Standard_Integer surfit = 0; surfit < 2; surfit++) {
 
         Handle(GeomAdaptor_HSurface) aGASurface = (surfit == 0) ? theSurface1 : theSurface2;
@@ -3819,8 +3684,9 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
           }
           else {
             Standard_Real aPeriod     = (parit == 0) ? aGASurface->UPeriod() : aGASurface->VPeriod();
-            Standard_Real anoffset = 0.;
-            Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
+            Standard_Real anoffset, anAdjustPar;
+            GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary,
+                                           aPeriod, anAdjustPar, anoffset);
 
             bIsPointOnBoundary=
               IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, aResolution, bIsOnFirstBoundary);
@@ -3830,11 +3696,11 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
               nbboundaries++;
             }
             else {
-              //check neighbourhood of boundary
-              Standard_Real anEpsilon = aResolution * 100.;
-              Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
-              anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
-                
+            //check neighbourhood of boundary
+            Standard_Real anEpsilon = aResolution * 100.;
+            Standard_Real aPart = ( aupperboundary - alowerboundary ) * 0.1;
+            anEpsilon = ( anEpsilon > aPart ) ? aPart : anEpsilon;
+            
               bIsNearBoundary = IsPointOnBoundary(anAdjustPar, alowerboundary, aupperboundary, 
                                                   anEpsilon, bIsOnFirstBoundary);
 
@@ -3847,7 +3713,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
           gp_Pnt2d aPZone = (surfit == 0) ? aTanZoneS1->Value(zIt) : aTanZoneS2->Value(zIt);
           Standard_Real aZoneRadius = aTanZoneRadius->Value(zIt);
 
-          Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
+          Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
           const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
           Standard_Real nU1, nV1;
             
@@ -3883,8 +3749,9 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
             Standard_Real aupperboundary = (bIsUBoundary) ? umax : vmax;
             Standard_Real aPeriod     = (bIsUBoundary) ? aGASurface->UPeriod() : aGASurface->VPeriod();
             Standard_Real aParameter = (bIsUBoundary) ? U : V;
-            Standard_Real anoffset = 0.;
-            Standard_Real anAdjustPar = AdjustPeriodic(aParameter, alowerboundary, aupperboundary, aPeriod, anoffset);
+            Standard_Real anoffset, anAdjustPar;
+            GeomInt::AdjustPeriodic(aParameter, alowerboundary, aupperboundary, 
+                                           aPeriod, anAdjustPar, anoffset);
 
             Standard_Real adist = (bIsFirstBoundary) ? fabs(anAdjustPar - alowerboundary) : fabs(anAdjustPar - aupperboundary);
             Standard_Real anotherPar = (bIsFirstBoundary) ? (aupperboundary - adist) : (alowerboundary + adist);
@@ -3912,7 +3779,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
               bCheckAngle1 = Standard_True;
               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(anewU, anewV));
 
-              if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
+              if(aNewVec.SquareMagnitude() < gp::Resolution()) {
                 aNewP.SetValue((surfit == 0), anewU, anewV);
                 bCheckAngle1 = Standard_False;
               }
@@ -3921,7 +3788,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
               bCheckAngle2 = Standard_True;
               aNewVec = gp_Vec2d(gp_Pnt2d(nU1, nV1), gp_Pnt2d(U, V));
 
-              if(aNewVec.SquareMagnitude() < (gp::Resolution() * gp::Resolution())) {
+              if(aNewVec.SquareMagnitude() < gp::Resolution()) {
                 bCheckAngle2 = Standard_False;
               }
             }
@@ -3930,7 +3797,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
               // assume there are at least two points in line (see "if" above)
               Standard_Integer anindexother = aneighbourpointindex;
 
-              while((anindexother <= aListOfIndex.Last()) && (anindexother >= aListOfIndex.First())) {
+              while((anindexother <= iLast) && (anindexother >= iFirst)) {
                 anindexother = (j == 0) ? (anindexother + 1) : (anindexother - 1);
                 const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(anindexother);
                 Standard_Real nU2, nV2;
@@ -3941,7 +3808,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
                   aPrevNeighbourPoint.ParametersOnS2(nU2, nV2);
                 gp_Vec2d aVecOld(gp_Pnt2d(nU2, nV2), gp_Pnt2d(nU1, nV1));
 
-                if(aVecOld.SquareMagnitude() <= (gp::Resolution() * gp::Resolution())) {
+                if(aVecOld.SquareMagnitude() <= gp::Resolution()) {
                   continue;
                 }
                 else {
@@ -3994,7 +3861,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
             else
               anewpoint = gp_Pnt2d( u2, v2 );
             
-            Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
+            Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
             Standard_Real nU1, nV1;
             
@@ -4034,7 +3901,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
           }
           else {
 
-            Standard_Integer aneighbourpointindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
+            Standard_Integer aneighbourpointindex1 = (j == 0) ? iFirst : iLast;
             const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneighbourpointindex1);
             Standard_Real nU1, nV1;
 
@@ -4046,7 +3913,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
             gp_Pnt2d ap2(nU1, nV1);
             Standard_Integer aneighbourpointindex2 = aneighbourpointindex1;
 
-            while((aneighbourpointindex2 <= aListOfIndex.Last()) && (aneighbourpointindex2 >= aListOfIndex.First())) {
+            while((aneighbourpointindex2 <= iLast) && (aneighbourpointindex2 >= iFirst)) {
               aneighbourpointindex2 = (j == 0) ? (aneighbourpointindex2 + 1) : (aneighbourpointindex2 - 1);
               const IntSurf_PntOn2S& aPrevNeighbourPoint = theWLine->Point(aneighbourpointindex2);
               Standard_Real nU2, nV2;
@@ -4058,7 +3925,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
               ap2.SetX(nU2);
               ap2.SetY(nV2);
 
-              if(ap1.SquareDistance(ap2) > (gp::Resolution() * gp::Resolution())) {
+              if(ap1.SquareDistance(ap2) > gp::Resolution()) {
                 break;
               }
             }  
@@ -4084,7 +3951,7 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
 
                 //Correction of projected coordinates. Begin
                 //Note, it may be shifted on a period
-                Standard_Integer aneindex1 = (j == 0) ? aListOfIndex.First() : aListOfIndex.Last();
+                Standard_Integer aneindex1 = (j == 0) ? iFirst : iLast;
                 const IntSurf_PntOn2S& aNeighbourPoint = theWLine->Point(aneindex1);
                 Standard_Real nUn, nVn;
                 
@@ -4167,10 +4034,6 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
         continue;
       }
       const TColStd_ListOfInteger& aListOfIndex = anArrayOfLines.Value(i);
-
-      if(aListOfIndex.Extent() < 2) {
-        continue;
-      }
       const TColStd_ListOfInteger& aListOfFLIndex = anArrayOfLineEnds.Value(i);
       Standard_Boolean bhasfirstpoint = (aListOfFLIndex.Extent() == 2);
       Standard_Boolean bhaslastpoint = (aListOfFLIndex.Extent() == 2);
@@ -4182,25 +4045,31 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
       if(!bhaslastpoint && !aListOfFLIndex.IsEmpty()) {
         bhaslastpoint = (i != nblines);
       }
-      Standard_Boolean bIsFirstInside = ((ifprm >= aListOfIndex.First()) && (ifprm <= aListOfIndex.Last()));
-      Standard_Boolean bIsLastInside =  ((ilprm >= aListOfIndex.First()) && (ilprm <= aListOfIndex.Last()));
+      
+      Standard_Integer iFirst = aListOfIndex.First();
+      Standard_Integer iLast  = aListOfIndex.Last();
+      Standard_Boolean bIsFirstInside = ((ifprm >= iFirst) && (ifprm <= iLast));
+      Standard_Boolean bIsLastInside =  ((ilprm >= iFirst) && (ilprm <= iLast));
 
       if(!bIsFirstInside && !bIsLastInside) {
-        if((ifprm < aListOfIndex.First()) && (ilprm > aListOfIndex.Last())) {
+        if((ifprm < iFirst) && (ilprm > iLast)) {
           // append whole line, and boundaries if neccesary
           if(bhasfirstpoint) {
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
+            pit = aListOfFLIndex.First();
+            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
             aLineOn2S->Add(aP);
           }
           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
 
           for(; anIt.More(); anIt.Next()) {
-            const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
+            pit = anIt.Value();
+            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
             aLineOn2S->Add(aP);
           }
 
           if(bhaslastpoint) {
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
+            pit = aListOfFLIndex.Last();
+            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
             aLineOn2S->Add(aP);
           }
 
@@ -4235,9 +4104,10 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
         TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
 
         for(; anIt.More(); anIt.Next()) {
-          if((anIt.Value() < ifprm) || (anIt.Value() > ilprm))
+          pit = anIt.Value();
+          if((pit < ifprm) || (pit > ilprm))
             continue;
-          const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
+          const IntSurf_PntOn2S& aP = theWLine->Point(pit);
           aLineOn2S->Add(aP);
         }
       }
@@ -4248,14 +4118,16 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
 
           for(; anIt.More(); anIt.Next()) {
-            if(anIt.Value() < ifprm)
+            pit = anIt.Value();
+            if(pit < ifprm)
               continue;
-            const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
+            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
             aLineOn2S->Add(aP);
           }
 
           if(bhaslastpoint) {
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.Last());
+            pit = aListOfFLIndex.Last();
+            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
             aLineOn2S->Add(aP);
           }
           // check end of split line (end is almost always)
@@ -4285,15 +4157,17 @@ Standard_Boolean DecompositionOfWLine(const Handle(IntPatch_WLine)& theWLine,
         if(bIsLastInside) {
           // append points from first boundary point to ilprm
           if(bhasfirstpoint) {
-            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(aListOfFLIndex.First());
+            pit = aListOfFLIndex.First();
+            const IntSurf_PntOn2S& aP = aSeqOfPntOn2S->Value(pit);
             aLineOn2S->Add(aP);
           }
           TColStd_ListIteratorOfListOfInteger anIt(aListOfIndex);
 
           for(; anIt.More(); anIt.Next()) {
-            if(anIt.Value() > ilprm)
+            pit = anIt.Value();
+            if(pit > ilprm)
               continue;
-            const IntSurf_PntOn2S& aP = theWLine->Point(anIt.Value());
+            const IntSurf_PntOn2S& aP = theWLine->Point(pit);
             aLineOn2S->Add(aP);
           }
         }
@@ -4324,7 +4198,7 @@ Standard_Boolean ParameterOutOfBoundary(const Standard_Real       theParameter,
                                         const Standard_Real       theOtherParameter,
                                         const Standard_Boolean    bIncreasePar,
                                         Standard_Real&            theNewParameter,
-                                        const Handle(BOPInt_Context)& aContext)
+                                        const Handle(IntTools_Context)& aContext)
 {
   Standard_Boolean bIsComputed = Standard_False;
   theNewParameter = theParameter;
@@ -4440,9 +4314,38 @@ Standard_Boolean ApproxWithPCurves(const gp_Cylinder& theCyl,
 {
   Standard_Boolean bRes = Standard_True;
   Standard_Real R1 = theCyl.Radius(), R2 = theSph.Radius();
-
-  if(R1 < 2.*R2) return bRes;
-
+  //
+  {
+    Standard_Real aD2, aRc2, aEps;
+    gp_Pnt aApexSph;
+    //
+    aEps=1.E-7;
+    aRc2=R1*R1;
+    //
+    const gp_Ax3& aAx3Sph=theSph.Position();
+    const gp_Pnt& aLocSph=aAx3Sph.Location();
+    const gp_Dir& aDirSph=aAx3Sph.Direction();
+    //
+    const gp_Ax1& aAx1Cyl=theCyl.Axis();
+    gp_Lin aLinCyl(aAx1Cyl);
+    //
+    aApexSph.SetXYZ(aLocSph.XYZ()+R2*aDirSph.XYZ());
+    aD2=aLinCyl.SquareDistance(aApexSph);
+    if (fabs(aD2-aRc2)<aEps) {
+      return !bRes;
+    }
+    //
+    aApexSph.SetXYZ(aLocSph.XYZ()-R2*aDirSph.XYZ());
+    aD2=aLinCyl.SquareDistance(aApexSph);
+    if (fabs(aD2-aRc2)<aEps) {
+      return !bRes;
+    }
+  }
+  //
+    
+  if(R1 < 2.*R2) {
+    return bRes;
+  }
   gp_Lin anCylAx(theCyl.Axis());
 
   Standard_Real aDist = anCylAx.Distance(theSph.Location());
@@ -4858,7 +4761,7 @@ Standard_Integer IndexType(const GeomAbs_SurfaceType aType)
   } 
   return aIndex;
 }
-#ifdef DEB_DUMPWLINE
+#ifdef OCCT_DEBUG_DUMPWLINE
 //=======================================================================
 //function : DumpWLine
 //purpose  : 
@@ -4913,127 +4816,31 @@ void RefineVector(gp_Vec2d& aV2D)
   }
   aV2D.SetCoord(aC[0], aC[1]);
 }
+
 //=======================================================================
-//function : FindMaxSquareDistance
-//purpose  : 
+// Function : MaxDistance
+// 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 MaxDistance(const Handle(Geom_Curve)& theC,
+                          const Standard_Real aT,
+                          GeomAPI_ProjectPointOnSurf& theProjPS)
 {
-  Standard_Real aA, aB, aCf, aX1, aX2, aF1, aF2, aX, aF;
+  Standard_Real aD;
+  gp_Pnt aP;
   //
-  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);
+  theC->D0(aT, aP);
+  theProjPS.Perform(aP);
+  aD = theProjPS.NbPoints() ? theProjPS.LowerDistance() : 0.;
   //
-  for(;;) {
-    //
-    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;
+  return aD;
 }
 
 //=======================================================================
 //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) 
+  Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC, 
+                               const TopoDS_Face& aFace) 
 {
   const Standard_Integer NPoints = 23;
   Standard_Integer i;
@@ -5044,7 +4851,7 @@ Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
   Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
   Standard_Real fp = aPC->FirstParameter();
   Standard_Real lp = aPC->LastParameter();
-  
+
 
   // adjust domain for periodic surfaces
   TopLoc_Location aLoc;
@@ -5093,11 +4900,37 @@ Standard_Boolean CheckPCurve(const Handle(Geom2d_Curve)& aPC,
       aT=aT+dT;
       aGAC.D0(aT, aP2D);
       aP2D.Coord(u,v);
-      if (umin-u > tolU || u-umax > tolU ||
+    if (umin-u > tolU || u-umax > tolU ||
           vmin-v > tolV || v-vmax > tolV) {
         return bRet;
-      }
-    }
+  }
+}
   }
   return !bRet;
 }
+//=======================================================================
+//function : CorrectPlaneBoundaries
+//purpose  : 
+//=======================================================================
+ void CorrectPlaneBoundaries(Standard_Real& aUmin,
+                             Standard_Real& aUmax, 
+                             Standard_Real& aVmin, 
+                             Standard_Real& aVmax) 
+{
+  if (!(Precision::IsInfinite(aUmin) ||
+        Precision::IsInfinite(aUmax))) { 
+    Standard_Real dU;
+    //
+    dU=0.1*(aUmax-aUmin);
+    aUmin=aUmin-dU;
+    aUmax=aUmax+dU;
+  }
+  if (!(Precision::IsInfinite(aVmin) ||
+        Precision::IsInfinite(aVmax))) { 
+    Standard_Real dV;
+    //
+    dV=0.1*(aVmax-aVmin);
+    aVmin=aVmin-dV;
+    aVmax=aVmax+dV;
+  }
+}