0023178: Intersection of cylinders fails to produce results
authornbv <nbv@opencascade.com>
Fri, 23 Sep 2016 14:24:52 +0000 (17:24 +0300)
committerapn <apn@opencascade.com>
Fri, 7 Oct 2016 10:37:31 +0000 (13:37 +0300)
1. Unification of trimmed and not-trimmed cylinders processing (IntPatch_Intersection::GeomGeomPerfomTrimSurf() method has been removed).
2. Interface of IntPatch_ImpImpIntersection::Perform(...) method has been changed.
3. Now, WLine purging is forbidden for Geom-Geom-Intersection.
4. Bnd_Range class has been created. See Bnd_Range.hxx for detail information.
5. Algorithm of AddBoundaryPoint function has been improved in order to obtain intersection points in both boundaries (VFirst and VLast of every surface).
6. Earlier, method Geom2dConvert::ConcatG1(...) increased resulted B-spline degree (in case of not succession of previous iteration). Now increased value has been limited by Geom2d_BSplineCurve::MaxDegree() value (max degree = 25).
7. Algorithm of B-spline closure definition has been changed in the methods Geom2dConvert::C0BSplineToC1BSplineCurve(...) and Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(...).

Creation of test case for this issue.
Adjusting test cases according to their new behavior.

Small correction in the code according to KGV's remark.

18 files changed:
src/BSplCLib/BSplCLib_2.cxx
src/Bnd/Bnd_Range.cxx [new file with mode: 0644]
src/Bnd/Bnd_Range.hxx [new file with mode: 0644]
src/Bnd/FILES
src/Geom2dConvert/Geom2dConvert.cxx
src/IntPatch/IntPatch_ImpImpIntersection.hxx
src/IntPatch/IntPatch_ImpImpIntersection_1.gxx
src/IntPatch/IntPatch_ImpImpIntersection_2.gxx
src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
src/IntPatch/IntPatch_Intersection.cxx
src/IntPatch/IntPatch_Intersection.hxx
src/IntPatch/IntPatch_WLineTool.cxx
src/IntPatch/IntPatch_WLineTool.hxx
tests/blend/simple/Q6
tests/bugs/modalg_5/bug25742_2
tests/bugs/modalg_6/bug23178 [new file with mode: 0644]
tests/bugs/modalg_6/bug27310_2
tests/bugs/modalg_6/bug27766

index 41b0835..6d7cbad 100644 (file)
@@ -31,8 +31,6 @@
 
 #include <math_Matrix.hxx>
 
-static const Standard_Integer aGlobalMaxDegree = 25;
-
 //=======================================================================
 //struct : BSplCLib_DataContainer 
 //purpose: Auxiliary structure providing buffers for poles and knots used in
@@ -44,8 +42,7 @@ struct BSplCLib_DataContainer
   BSplCLib_DataContainer(Standard_Integer Degree)
   {
     (void)Degree; // avoid compiler warning
-    Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
-        BSplCLib::MaxDegree() > aGlobalMaxDegree,
+    Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree(),
         "BSplCLib: bspline degree is greater than maximum supported");
   }
 
@@ -326,7 +323,7 @@ BSplCLib::BuildBSpMatrix(const  TColStd_Array1OfReal&     Parameters,
   ReturnCode = 0,
   FirstNonZeroBsplineIndex,
   BandWidth,
-  MaxOrder = aGlobalMaxDegree+1,
+  MaxOrder = BSplCLib::MaxDegree() + 1,
   Order ;
   
   math_Matrix   BSplineBasis(1, MaxOrder,
diff --git a/src/Bnd/Bnd_Range.cxx b/src/Bnd/Bnd_Range.cxx
new file mode 100644 (file)
index 0000000..d79b211
--- /dev/null
@@ -0,0 +1,37 @@
+// Created on: 2016-06-07
+// Created by: Nikolai BUKHALOV
+// Copyright (c) 2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+
+#include <Bnd_Range.hxx>
+
+//=======================================================================
+//function : Common
+//purpose  : 
+//=======================================================================
+void Bnd_Range::Common(const Bnd_Range& theOther)
+{
+  if(theOther.IsVoid())
+  {
+    SetVoid();
+  }
+
+  if(IsVoid())
+  {
+    return;
+  }
+
+  myFirst = Max(myFirst, theOther.myFirst);
+  myLast = Min(myLast, theOther.myLast);
+}
diff --git a/src/Bnd/Bnd_Range.hxx b/src/Bnd/Bnd_Range.hxx
new file mode 100644 (file)
index 0000000..439008f
--- /dev/null
@@ -0,0 +1,124 @@
+// Created on: 2016-06-07
+// Created by: Nikolai BUKHALOV
+// Copyright (c) 2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _Bnd_Range_HeaderFile
+#define _Bnd_Range_HeaderFile
+
+#include <Standard_Real.hxx>
+#include <Standard_ConstructionError.hxx>
+
+//! This class describes a range in 1D space restricted
+//! by two real values.
+//! A range can be void indicating there is no point included in the range.
+
+class Bnd_Range
+{
+public:
+
+  //! Default constructor. Creates VOID range.
+  Bnd_Range() : myFirst(0.0), myLast(-1.0)
+  {
+  };
+
+  //! Constructor. Never creates VOID range.
+  Bnd_Range(const Standard_Real theMin, const Standard_Real theMax) : 
+                                                    myFirst(theMin), myLast(theMax)
+  {
+    if(myLast < myFirst)
+      Standard_ConstructionError::Raise("Last < First");
+  };
+
+  //! Replaces <this> with common-part of <this> and theOther
+  Standard_EXPORT void Common(const Bnd_Range& theOther);
+
+  //! Extends <this> to include theParameter
+  void Add(const Standard_Real theParameter)
+  {
+    if(IsVoid())
+    {
+      myFirst = myLast = theParameter;
+      return;
+    }
+
+    myFirst = Min(myFirst, theParameter);
+    myLast = Max(myLast, theParameter);
+  }
+
+  //! Obtain MIN boundary of <this>.
+  //! If <this> is VOID the method returns false.
+  Standard_Boolean GetMin(Standard_Real& thePar) const
+  {
+    if(IsVoid())
+    {
+      return Standard_False;
+    }
+
+    thePar = myFirst;
+    return Standard_True;
+  }
+
+  //! Obtain MAX boundary of <this>.
+  //! If <this> is VOID the method returns false.
+  Standard_Boolean GetMAX(Standard_Real& thePar) const
+  {
+    if(IsVoid())
+    {
+      return Standard_False;
+    }
+
+    thePar = myLast;
+    return Standard_True;
+  }
+
+  //! Returns range value (MAX-MIN). Returns negative value for VOID range.
+  Standard_Real Delta() const
+  {
+    return (myLast - myFirst);
+  }
+
+  //! Is <this> initialized.
+  Standard_Boolean IsVoid() const
+  {
+    return (myLast < myFirst);
+  }
+
+  //! Initializes <this> by default parameters. Makes <this> VOID.
+  void SetVoid()
+  {
+    myLast = -1.0;
+    myFirst = 0.0;
+  }
+
+  //! Extends this to the given value (in both side)
+  void Enlarge(const Standard_Real theDelta)
+  {
+    if (IsVoid())
+    {
+      return;
+    }
+
+    myFirst -= theDelta;
+    myLast += theDelta;
+  }
+
+private:
+  //! Start of range
+  Standard_Real myFirst;
+
+  //! End of range
+  Standard_Real myLast;
+};
+
+#endif
\ No newline at end of file
index 6c20deb..845b525 100644 (file)
@@ -24,6 +24,8 @@ Bnd_Box2d.hxx
 Bnd_HArray1OfBox.hxx
 Bnd_HArray1OfBox2d.hxx
 Bnd_HArray1OfSphere.hxx
+Bnd_Range.cxx
+Bnd_Range.hxx
 Bnd_SeqOfBox.hxx
 Bnd_Sphere.cxx
 Bnd_Sphere.hxx
index 8b952ac..9190434 100644 (file)
@@ -1025,6 +1025,7 @@ void  Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
    for (j=1;j<=nb_curve-1;j++){                //boucle secondaire a l'interieur de chaque groupe
      Curve1=ArrayOfCurves(j);
      if ( (j==(nb_curve-1)) &&(Need2DegRepara(ArrayOfCurves))){ 
+       const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree());
        Curve2->D1(Curve2->LastParameter(),Pint,Vec1);
        Curve1->D1(Curve1->FirstParameter(),Pint,Vec2);
        lambda=Vec2.Magnitude()/Vec1.Magnitude();
@@ -1074,7 +1075,7 @@ void  Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
                                        Curve1FlatKnots,
                                        Curve1Poles,
                                        FlatKnots,
-                                       2*Curve1->Degree(),
+                                        aNewCurveDegree,
                                        NewPoles,
                                        Status
                                        );
@@ -1084,7 +1085,7 @@ void  Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
                                        Curve1FlatKnots,
                                        Curve1Weights,
                                        FlatKnots,
-                                       2*Curve1->Degree(),
+                                        aNewCurveDegree,
                                        NewWeights,
                                        Status
                                        );
@@ -1110,7 +1111,7 @@ void  Geom2dConvert::ConcatG1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
        for (ii=1;ii<=NewPoles.Length();ii++)
         for (jj=1;jj<=2;jj++)
           NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii));
-       Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree());
+       Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree);
      }
      Geom2dConvert_CompCurveToBSplineCurve C(Curve2);
      fusion=C.Add(Curve1,
@@ -1276,6 +1277,8 @@ void  Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
      else
        Curve1=ArrayOfCurves(j);
      
+     const Standard_Integer aNewCurveDegree = Min(2 * Curve1->Degree(), Geom2d_BSplineCurve::MaxDegree());
+
      if (j==0)                                      //initialisation en debut de groupe
        Curve2=Curve1;
      else{
@@ -1315,7 +1318,7 @@ void  Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
         TColStd_Array1OfReal FlatKnots(1,Curve1FlatKnots.Length()+(Curve1->Degree()*Curve1->NbKnots()));
         
         BSplCLib::KnotSequence(KnotC1,KnotC1Mults,FlatKnots);
-        TColgp_Array1OfPnt2d  NewPoles(1,FlatKnots.Length()-(2*Curve1->Degree()+1));
+         TColgp_Array1OfPnt2d  NewPoles(1, FlatKnots.Length() - (aNewCurveDegree + 1));
         Standard_Integer      Status;
         TColStd_Array1OfReal Curve1Weights(1,Curve1->NbPoles());
         Curve1->Weights(Curve1Weights);
@@ -1330,18 +1333,18 @@ void  Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
                                          Curve1FlatKnots,
                                          Curve1Poles,
                                          FlatKnots,
-                                         2*Curve1->Degree(),
+                                          aNewCurveDegree,
                                          NewPoles,
                                          Status
                                          );
-        TColStd_Array1OfReal NewWeights(1,FlatKnots.Length()-(2*Curve1->Degree()+1));
+         TColStd_Array1OfReal NewWeights(1, FlatKnots.Length() - (aNewCurveDegree + 1));
 //      BSplCLib::FunctionReparameterise(reparameterise_evaluator,
         BSplCLib::FunctionReparameterise(ev,
                                          Curve1->Degree(),
                                          Curve1FlatKnots,
                                          Curve1Weights,
                                          FlatKnots,
-                                         2*Curve1->Degree(),
+                                          aNewCurveDegree,
                                          NewWeights,
                                          Status
                                          );
@@ -1349,7 +1352,7 @@ void  Geom2dConvert::ConcatC1(TColGeom2d_Array1OfBSplineCurve&           ArrayOf
           for (jj=1;jj<=2;jj++)
             NewPoles(ii).SetCoord(jj,NewPoles(ii).Coord(jj)/NewWeights(ii));
         }
-        Curve1= new Geom2d_BSplineCurve(NewPoles,NewWeights,KnotC1,KnotC1Mults,2*Curve1->Degree());
+         Curve1 = new Geom2d_BSplineCurve(NewPoles, NewWeights, KnotC1, KnotC1Mults, aNewCurveDegree);
        }
        Geom2dConvert_CompCurveToBSplineCurve C(Curve2);
        fusion=C.Add(Curve1,
@@ -1420,7 +1423,7 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS,
  Standard_Integer                 i,j,nbcurveC1=1;
  Standard_Real                    U1,U2;
  Standard_Boolean                 closed_flag = Standard_False ;
- gp_Pnt2d                         point;
+ gp_Pnt2d                         point1, point2;
  gp_Vec2d                         V1,V2;
  Standard_Boolean                 fusion;
 
@@ -1453,15 +1456,20 @@ void Geom2dConvert::C0BSplineToC1BSplineCurve(Handle(Geom2d_BSplineCurve)& BS,
      BSbis->Segment(U1,U2);
      ArrayOfCurves(i)=BSbis;
    }
+
+   const Standard_Real anAngularToler = 1.0e-7;
    Handle(TColStd_HArray1OfInteger) ArrayOfIndices;
    Handle(TColGeom2d_HArray1OfBSplineCurve) ArrayOfConcatenated;
     
-   BS->D1(BS->FirstParameter(),point,V1);  //a verifier
-   BS->D1(BS->LastParameter(),point,V2);
+   BS->D1(BS->FirstParameter(),point1,V1);  //a verifier
+   BS->D1(BS->LastParameter(),point2,V2);
    
-   if ((BS->IsClosed())&&(V1.IsParallel(V2,Precision::Confusion())))
-     closed_flag = Standard_True ;
-    
+   if ((point1.SquareDistance(point2) < gp::Resolution()) &&
+       (V1.IsParallel(V2, anAngularToler)))
+   {
+     closed_flag = Standard_True;
+   }
+
    Geom2dConvert::ConcatC1(ArrayOfCurves,
                           ArrayOfToler,
                           ArrayOfIndices,
@@ -1510,7 +1518,7 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline
   Standard_Integer                 i,j,nbcurveC1=1;
   Standard_Real                    U1,U2;
   Standard_Boolean                 closed_flag = Standard_False ;
-  gp_Pnt2d                         point;
+  gp_Pnt2d                         point1, point2;
   gp_Vec2d                         V1,V2;
 //  Standard_Boolean                 fusion;
   
@@ -1544,11 +1552,14 @@ void Geom2dConvert::C0BSplineToArrayOfC1BSplineCurve(const Handle(Geom2d_BSpline
     
     Handle(TColStd_HArray1OfInteger) ArrayOfIndices;
     
-    BS->D1(BS->FirstParameter(),point,V1);  
-    BS->D1(BS->LastParameter(),point,V2);
+    BS->D1(BS->FirstParameter(),point1,V1);  
+    BS->D1(BS->LastParameter(),point2,V2);
     
-    if ((BS->IsClosed())&&(V1.IsParallel(V2,AngularTolerance)))
-      closed_flag = Standard_True ;
+    if (((point1.SquareDistance(point2) < gp::Resolution())) &&
+        (V1.IsParallel(V2, AngularTolerance)))
+    {
+      closed_flag = Standard_True;
+    }
     
     Geom2dConvert::ConcatC1(ArrayOfCurves,
                            ArrayOfToler,
index 23153eb..ece87d6 100644 (file)
@@ -62,7 +62,14 @@ public:
   //! When intersection result returns IntPatch_RLine and another
   //! IntPatch_Line (not restriction) we (in case of theIsReqToKeepRLine==TRUE)
   //! will always keep both lines even if they are coincided.
-  Standard_EXPORT void Perform (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, const Standard_Boolean isTheTrimmed = Standard_False, const Standard_Boolean theIsReqToKeepRLine = Standard_False);
+  Standard_EXPORT void Perform (const Handle(Adaptor3d_HSurface)& S1,
+                                const Handle(Adaptor3d_TopolTool)& D1,
+                                const Handle(Adaptor3d_HSurface)& S2,
+                                const Handle(Adaptor3d_TopolTool)& D2,
+                                const Standard_Real TolArc,
+                                const Standard_Real TolTang,
+                                const Standard_Boolean theIsReqToKeepRLine =
+                                                                  Standard_False);
   
   //! Returns True if the calculus was succesfull.
     Standard_Boolean IsDone() const;
@@ -71,8 +78,8 @@ public:
     Standard_Boolean IsEmpty() const;
   
   //! Returns True if the two patches are considered as
-  //! entierly tangent, i-e every restriction arc of one
-  //! patch is inside the geometric base of the otehr patch.
+  //! entirely tangent, i.e every restriction arc of one
+  //! patch is inside the geometric base of the other patch.
     Standard_Boolean TangentFaces() const;
   
   //! Returns True when the TangentFaces returns True and the
index 18146a1..e2f049e 100644 (file)
@@ -70,16 +70,7 @@ static void ProcessBounds(const Handle(IntPatch_ALine)&,
                          const Standard_Real);
 
 
-static Standard_Boolean IntCyCy(const IntSurf_Quadric&,
-                               const IntSurf_Quadric&,
-                               const Standard_Real,
-                               Standard_Boolean&,
-                               Standard_Boolean&,
-                               Standard_Boolean&,
-                               IntPatch_SequenceOfLine&,
-                               IntPatch_SequenceOfPoint&);
-
-static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1,
+static Standard_Boolean IntCyCy(const IntSurf_Quadric& theQuad1,
                                     const IntSurf_Quadric& theQuad2,
                                     const Standard_Real theTol3D,
                                     const Standard_Real theTol2D,
@@ -87,6 +78,8 @@ static Standard_Boolean IntCyCyTrim(const IntSurf_Quadric& theQuad1,
                                     const Bnd_Box2d& theUVSurf2,
                                     const Standard_Boolean isTheReverse,
                                     Standard_Boolean& isTheEmpty,
+                                    Standard_Boolean& isTheSameSurface,
+                                    Standard_Boolean& isTheMultiplePoint,
                                     IntPatch_SequenceOfLine& theSlin,
                                     IntPatch_SequenceOfPoint& theSPnt);
 
index b1d3dfc..5972adf 100644 (file)
@@ -14,6 +14,8 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+#include <IntPatch_WLine.hxx>
+
 static 
   Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS,
                            GeomAbs_SurfaceType& theTS,
@@ -40,7 +42,7 @@ IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection
         const Standard_Real TolTang,
         const Standard_Boolean theIsReqToKeepRLine)
 {
-  Perform(S1,D1,S2,D2,TolArc,TolTang, Standard_False, theIsReqToKeepRLine);
+  Perform(S1,D1,S2,D2,TolArc,TolTang, theIsReqToKeepRLine);
 }
 //=======================================================================
 //function : Perform
@@ -52,13 +54,14 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
                                           const Handle(Adaptor3d_TopolTool)& D2,
                                           const Standard_Real TolArc,
                                           const Standard_Real TolTang,
-                                          const Standard_Boolean isTheTrimmed,
-                                          const Standard_Boolean theIsReqToKeepRLine) {
+                                          const Standard_Boolean theIsReqToKeepRLine)
+{
   done = Standard_False;
-  Standard_Boolean isTrimmed = isTheTrimmed;
   spnt.Clear();
   slin.Clear();
 
+  Standard_Boolean isPostProcessingRequired = Standard_True;
+
   empt = Standard_True;
   tgte = Standard_False;
   oppo = Standard_False;
@@ -149,59 +152,53 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
     case 22:
       { // Cylinder/Cylinder
         Standard_Boolean isDONE = Standard_False;
-        
-        if(!isTrimmed)
+        Bnd_Box2d aBox1, aBox2;
+
+        const Standard_Real aU1f = S1->FirstUParameter();
+        Standard_Real aU1l = S1->LastUParameter();
+        const Standard_Real aU2f = S2->FirstUParameter();
+        Standard_Real aU2l = S2->LastUParameter();
+
+        const Standard_Real anUperiod = 2.0*M_PI;
+
+        if(aU1l - aU1f > anUperiod)
+          aU1l = aU1f + anUperiod;
+
+        if(aU2l - aU2f > anUperiod)
+          aU2l = aU2f + anUperiod;
+
+        aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter()));
+        aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter()));
+        aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter()));
+        aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter()));
+
+        // Resolution is too big if the cylinder radius is
+        // too small. Therefore, we shall bounde its value above. 
+        // Here, we use simple constant.
+        const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang),
+                                            S2->UResolution(TolTang)));
+
+        //The bigger range the bigger nunber of points in Walking-line (WLine)
+        //we will be able to add and, consequently, we will obtain more
+        //precise intersection line.
+        //Every point of WLine is determined as function from U1-parameter,
+        //where U1 is U-parameter on 1st quadric.
+        //Therefore, we should use quadric with bigger range as 1st parameter
+        //in IntCyCy() function.
+        //On the other hand, there is no point in reversing in case of
+        //analytical intersection (when result is line, ellipse, point...).
+        //This result is independent of the arguments order.
+        Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f));
+
+        if(isReversed)
         {
-          isDONE = IntCyCy(quad1, quad2, TolTang, empt,
-                            SameSurf, multpoint, slin, spnt);
+          isDONE = IntCyCy(quad2, quad1, TolTang, a2DTol, aBox2, aBox1,
+                                    Standard_True, empt, SameSurf, multpoint, slin, spnt);
         }
         else
         {
-          Bnd_Box2d aBox1, aBox2;
-
-          const Standard_Real aU1f = S1->FirstUParameter();
-          Standard_Real aU1l = S1->LastUParameter();
-          const Standard_Real aU2f = S2->FirstUParameter();
-          Standard_Real aU2l = S2->LastUParameter();
-
-          const Standard_Real anUperiod = 2.0*M_PI;
-
-          if(aU1l - aU1f > anUperiod)
-            aU1l = aU1f + anUperiod;
-
-          if(aU2l - aU2f > anUperiod)
-            aU2l = aU2f + anUperiod;
-
-          aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter()));
-          aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter()));
-          aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter()));
-          aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter()));
-
-          // Resolution is too big if the cylinder radius is
-          // too small. Therefore, we shall bounde its value above. 
-          // Here, we use simple constant.
-          const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang),
-                                              S2->UResolution(TolTang)));
-
-          Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f));
-
-          if(isReversed)
-          {
-            isDONE = IntCyCyTrim(quad2, quad1, TolTang, a2DTol, aBox2, aBox1,
-                                                  isReversed, empt, slin, spnt);
-          }
-          else
-          {
-            isDONE = IntCyCyTrim(quad1, quad2, TolTang, a2DTol, aBox1, aBox2,
-                                                  isReversed, empt, slin, spnt);
-          }
-
-          if(!isDONE)
-          {
-            isDONE = IntCyCy(quad1, quad2, TolTang, empt,
-                            SameSurf, multpoint, slin, spnt);
-            isTrimmed = Standard_False;
-          }
+          isDONE = IntCyCy(quad1, quad2, TolTang, a2DTol, aBox1, aBox2,
+                                    Standard_False, empt, SameSurf, multpoint, slin, spnt);
         }
 
         if (!isDONE)
@@ -210,6 +207,17 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
         }
 
         bEmpty = empt;
+        if(!slin.IsEmpty())
+        {
+          const Handle(IntPatch_WLine)& aWLine = 
+                                    Handle(IntPatch_WLine)::DownCast(slin.Value(1));
+
+          if(!aWLine.IsNull())
+          {//No geometric solution
+            isPostProcessingRequired = Standard_False;
+          }
+        }
+
         break;
       }
     //
@@ -302,7 +310,7 @@ void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
   }
   //
 
-  if(!isTrimmed)
+  if(isPostProcessingRequired)
   {
     if (!SameSurf) {
       AFunc.SetQuadric(quad2);
index cee3358..f5aca30 100644 (file)
 // commercial license or contractual agreement.
 
 #include <algorithm>
-#include <Standard_DivideByZero.hxx>
-#include <IntAna_ListOfCurve.hxx>
+#include <Bnd_Range.hxx>
 #include <IntAna_ListIteratorOfListOfCurve.hxx>
+#include <IntAna_ListOfCurve.hxx>
 #include <IntPatch_WLine.hxx>
-
+#include <Standard_DivideByZero.hxx>
 #include <math_Matrix.hxx>
 
 //If Abs(a) <= aNulValue then it is considered that a = 0.
 static const Standard_Real aNulValue = 1.0e-11;
 
-struct stCoeffsValue;
+static void ShortCosForm( const Standard_Real theCosFactor,
+                          const Standard_Real theSinFactor,
+                          Standard_Real& theCoeff,
+                          Standard_Real& theAngle);
 //
 static 
   Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
@@ -33,10 +36,6 @@ static
                                IntAna_Curve& aC,
                                const Standard_Real aTol,
                                IntAna_ListOfCurve& aLC);
-static
-  Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
-                              const gp_Cylinder& Cy2,
-                              const Standard_Real Tol);
 
 static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
                                       const Standard_Real theUlTarget,
@@ -45,16 +44,399 @@ static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
                                       const Standard_Real thePeriod,
                                       const Standard_Boolean theFlForce);
 
+
+class ComputationMethods
+{
+public:
+  //Stores equations coefficients
+  struct stCoeffsValue
+  {
+    stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
+
+    math_Vector mVecA1;
+    math_Vector mVecA2;
+    math_Vector mVecB1;
+    math_Vector mVecB2;
+    math_Vector mVecC1;
+    math_Vector mVecC2;
+    math_Vector mVecD;
+
+    Standard_Real mK21; //sinU2
+    Standard_Real mK11; //sinU1
+    Standard_Real mL21; //cosU2
+    Standard_Real mL11;  //cosU1
+    Standard_Real mM1;  //Free member
+
+    Standard_Real mK22; //sinU2
+    Standard_Real mK12; //sinU1
+    Standard_Real mL22; //cosU2
+    Standard_Real mL12; //cosU1
+    Standard_Real mM2; //Free member
+
+    Standard_Real mK1;
+    Standard_Real mL1;
+    Standard_Real mK2;
+    Standard_Real mL2;
+
+    Standard_Real mFIV1;
+    Standard_Real mPSIV1;
+    Standard_Real mFIV2;
+    Standard_Real mPSIV2;
+
+    Standard_Real mB;
+    Standard_Real mC;
+    Standard_Real mFI1;
+    Standard_Real mFI2;
+  };
+
+
+  //! Determines, if U2(U1) function is increasing.
+  static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
+                                             const Standard_Integer theWLIndex,
+                                             const stCoeffsValue& theCoeffs,
+                                             const Standard_Real thePeriod,
+                                             Standard_Boolean& theIsIncreasing);
+
+  //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
+  //! esimates the tolerance of U2-computing (estimation result is
+  //! assigned to *theDelta value).
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+                                                const Standard_Integer theWLIndex,
+                                                const stCoeffsValue& theCoeffs,
+                                                Standard_Real& theU2,
+                                                Standard_Real* const theDelta = 0);
+
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
+                                                  const Standard_Real theU2,
+                                                  const stCoeffsValue& theCoeffs,
+                                                  Standard_Real& theV1,
+                                                  Standard_Real& theV2);
+
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+                                                  const Standard_Integer theWLIndex,
+                                                  const stCoeffsValue& theCoeffs,
+                                                  Standard_Real& theU2,
+                                                  Standard_Real& theV1,
+                                                  Standard_Real& theV2);
+
+};
+
+ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
+                                                 const gp_Cylinder& theCyl2):
+  mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
+  mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
+  mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
+  mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
+  mVecC1(theCyl1.Axis().Direction().XYZ()),
+  mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
+  mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
+{
+  enum CoupleOfEquation
+  {
+    COENONE = 0,
+    COE12 = 1,
+    COE23 = 2,
+    COE13 = 3
+  }aFoundCouple = COENONE;
+
+
+  Standard_Real aDetV1V2 = 0.0;
+
+  const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
+  const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
+  const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
+  const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
+  const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
+  const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
+
+  if(anAbsD1 >= anAbsD2)
+  {
+    if(anAbsD3 > anAbsD1)
+    {
+      aFoundCouple = COE13;
+      aDetV1V2 = aDelta3;
+    }
+    else
+    {
+      aFoundCouple = COE12;
+      aDetV1V2 = aDelta1;
+    }
+  }
+  else
+  {
+    if(anAbsD3 > anAbsD2)
+    {
+      aFoundCouple = COE13;
+      aDetV1V2 = aDelta3;
+    }
+    else
+    {
+      aFoundCouple = COE23;
+      aDetV1V2 = aDelta2;
+    }
+  }
+
+  // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
+  // cross-product between directions (i.e. sine of angle).
+  // If sine is too small then sine is (approx.) equal to angle itself.
+  // Therefore, in this case we should compare sine with angular tolerance. 
+  // This constant is used for check if axes are parallel (see constructor
+  // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
+  if(Abs(aDetV1V2) < Precision::Angular())
+  {
+    Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
+  }
+
+  switch(aFoundCouple)
+  {
+  case COE12:
+    break;
+  case COE23:
+    {
+      math_Vector aVTemp(mVecA1);
+      mVecA1(1) = aVTemp(2);
+      mVecA1(2) = aVTemp(3);
+      mVecA1(3) = aVTemp(1);
+
+      aVTemp = mVecA2;
+      mVecA2(1) = aVTemp(2);
+      mVecA2(2) = aVTemp(3);
+      mVecA2(3) = aVTemp(1);
+
+      aVTemp = mVecB1;
+      mVecB1(1) = aVTemp(2);
+      mVecB1(2) = aVTemp(3);
+      mVecB1(3) = aVTemp(1);
+
+      aVTemp = mVecB2;
+      mVecB2(1) = aVTemp(2);
+      mVecB2(2) = aVTemp(3);
+      mVecB2(3) = aVTemp(1);
+
+      aVTemp = mVecC1;
+      mVecC1(1) = aVTemp(2);
+      mVecC1(2) = aVTemp(3);
+      mVecC1(3) = aVTemp(1);
+
+      aVTemp = mVecC2;
+      mVecC2(1) = aVTemp(2);
+      mVecC2(2) = aVTemp(3);
+      mVecC2(3) = aVTemp(1);
+
+      aVTemp = mVecD;
+      mVecD(1) = aVTemp(2);
+      mVecD(2) = aVTemp(3);
+      mVecD(3) = aVTemp(1);
+
+      break;
+    }
+  case COE13:
+    {
+      math_Vector aVTemp = mVecA1;
+      mVecA1(2) = aVTemp(3);
+      mVecA1(3) = aVTemp(2);
+
+      aVTemp = mVecA2;
+      mVecA2(2) = aVTemp(3);
+      mVecA2(3) = aVTemp(2);
+
+      aVTemp = mVecB1;
+      mVecB1(2) = aVTemp(3);
+      mVecB1(3) = aVTemp(2);
+
+      aVTemp = mVecB2;
+      mVecB2(2) = aVTemp(3);
+      mVecB2(3) = aVTemp(2);
+
+      aVTemp = mVecC1;
+      mVecC1(2) = aVTemp(3);
+      mVecC1(3) = aVTemp(2);
+
+      aVTemp = mVecC2;
+      mVecC2(2) = aVTemp(3);
+      mVecC2(3) = aVTemp(2);
+
+      aVTemp = mVecD;
+      mVecD(2) = aVTemp(3);
+      mVecD(3) = aVTemp(2);
+
+      break;
+    }
+  default:
+    break;
+  }
+
+  //------- For V1 (begin)
+  //sinU2
+  mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
+  //sinU1
+  mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
+  //cosU2
+  mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
+  //cosU1
+  mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
+  //Free member
+  mM1 =  (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
+  //------- For V1 (end)
+
+  //------- For V2 (begin)
+  //sinU2
+  mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
+  //sinU1
+  mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
+  //cosU2
+  mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
+  //cosU1
+  mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
+  //Free member
+  mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
+  //------- For V1 (end)
+
+  ShortCosForm(mL11, mK11, mK1, mFIV1);
+  ShortCosForm(mL21, mK21, mL1, mPSIV1);
+  ShortCosForm(mL12, mK12, mK2, mFIV2);
+  ShortCosForm(mL22, mK22, mL2, mPSIV2);
+
+  const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
+    aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
+    aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
+    aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
+
+  mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2;  //Free
+
+  Standard_Real aA = 0.0;
+
+  ShortCosForm(aB2,aB1,mB,mFI1);
+  ShortCosForm(aA2,aA1,aA,mFI2);
+
+  mB /= aA;
+  mC /= aA;
+}
+
+class WorkWithBoundaries
+{
+public:
+  enum SearchBoundType
+  {
+    SearchNONE = 0,
+    SearchV1 = 1,
+    SearchV2 = 2
+  };
+
+  struct StPInfo
+  {
+    StPInfo()
+    {
+      mySurfID = 0;
+      myU1 = RealLast();
+      myV1 = RealLast();
+      myU2 = RealLast();
+      myV2 = RealLast();
+    }
+
+    //Equal to 0 for 1st surface non-zerro for 2nd one.
+    Standard_Integer mySurfID;
+
+    Standard_Real myU1;
+    Standard_Real myV1;
+    Standard_Real myU2;
+    Standard_Real myV2;
+
+    bool operator>(const StPInfo& theOther) const
+    {
+      return myU1 > theOther.myU1;
+    }
+
+    bool operator<(const StPInfo& theOther) const
+    {
+      return myU1 < theOther.myU1;
+    }
+
+    bool operator==(const StPInfo& theOther) const
+    {
+      return myU1 == theOther.myU1;
+    }
+  };
+
+  WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
+                     const IntSurf_Quadric& theQuad2,
+                     const ComputationMethods::stCoeffsValue& theCoeffs,
+                     const Bnd_Box2d& theUVSurf1,
+                     const Bnd_Box2d& theUVSurf2,
+                     const Standard_Integer theNbWLines,
+                     const Standard_Real thePeriod,
+                     const Standard_Real theTol3D,
+                     const Standard_Real theTol2D,
+                     const Standard_Boolean isTheReverse) :
+      myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
+      myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
+      myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
+      myIsReverse(isTheReverse)
+  {
+  };
+
+  void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
+                        const Standard_Real theU1,
+                        const Standard_Real theU2,
+                        const Standard_Real theV1,
+                        const Standard_Real theV1Prev,
+                        const Standard_Real theV2,
+                        const Standard_Real theV2Prev,
+                        const Standard_Integer theWLIndex,
+                        const Standard_Boolean theFlForce,
+                        Standard_Boolean& isTheFound1,
+                        Standard_Boolean& isTheFound2) const;
+  
+  Standard_Boolean BoundariesComputing(Standard_Real theU1f[],
+                                       Standard_Real theU1l[]) const;
+  
+  void BoundaryEstimation(const gp_Cylinder& theCy1,
+                          const gp_Cylinder& theCy2,
+                          Bnd_Range& theOutBoxS1,
+                          Bnd_Range& theOutBoxS2) const;
+
+protected:
+
+  Standard_Boolean
+              SearchOnVBounds(const SearchBoundType theSBType,
+                              const Standard_Real theVzad,
+                              const Standard_Real theVInit,
+                              const Standard_Real theInitU2,
+                              const Standard_Real theInitMainVar,
+                              Standard_Real& theMainVariableValue) const;
+
+  const WorkWithBoundaries& operator=(const WorkWithBoundaries&);
+
+private:
+  friend class ComputationMethods;
+
+  const IntSurf_Quadric& myQuad1;
+  const IntSurf_Quadric& myQuad2;
+  const ComputationMethods::stCoeffsValue& myCoeffs;
+  const Bnd_Box2d& myUVSurf1;
+  const Bnd_Box2d& myUVSurf2;
+  const Standard_Integer myNbWLines;
+  const Standard_Real myPeriod;
+  const Standard_Real myTol3D;
+  const Standard_Real myTol2D;
+  const Standard_Boolean myIsReverse;
+};
+
+static 
+  Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
+                               const gp_Cone& aCo,
+                               IntAna_Curve& aC,
+                               const Standard_Real aTol,
+                               IntAna_ListOfCurve& aLC);
+
 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const IntSurf_Quadric& theQuad2,
                                   const Handle(IntSurf_LineOn2S)& theLine,
-                                  const stCoeffsValue& theCoeffs,
+                                  const ComputationMethods::stCoeffsValue& theCoeffs,
                                   const Standard_Integer theWLIndex,
                                   const Standard_Integer theMinNbPoints,
                                   const Standard_Integer theStartPointOnLine,
                                   const Standard_Integer theEndPointOnLine,
-                                  const Standard_Real theU2f,
-                                  const Standard_Real theU2l,
                                   const Standard_Real theTol2D,
                                   const Standard_Real thePeriodOfSurf2,
                                   const Standard_Boolean isTheReverse);
@@ -75,6 +457,29 @@ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
 }
 
 //=======================================================================
+//function : ExtremaLineLine
+//purpose  : Computes extrema between the given lines. Returns parameters
+//          on correspond curve. 
+//=======================================================================
+static inline void ExtremaLineLine(const gp_Ax1& theC1,
+                                   const gp_Ax1& theC2,
+                                   const Standard_Real theCosA,
+                                   const Standard_Real theSqSinA,
+                                   Standard_Real& thePar1,
+                                   Standard_Real& thePar2)
+{
+  const gp_Dir &aD1 = theC1.Direction(),
+               &aD2 = theC2.Direction();
+
+  const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
+  const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
+                      aD2L = aD2.XYZ().Dot(aL1L2);
+
+  thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
+  thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
+}
+
+//=======================================================================
 //function : VBoundaryPrecise
 //purpose  : By default, we shall consider, that V1 and V2 will increase
 //            if U1 increases. But if it is not, new V1set and/or V2set
@@ -165,7 +570,7 @@ static Standard_Boolean StepComputing(const math_Matrix& theMatr,
                                       Standard_Real& theV1Found,
                                       Standard_Real& theV2Found*/)
 {
-#ifdef OCCT_DEBUG
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
   bool flShow = false;
 
   if(flShow)
@@ -426,22 +831,24 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
     alig->SetLastPoint(alig->NbVertex());
   }
 }
+
 //=======================================================================
-//function : IntCyCy
+//function : CyCyAnalyticalIntersect
 //purpose  : 
 //=======================================================================
-Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
-                        const IntSurf_Quadric& Quad2,
-                        const Standard_Real Tol,
-                        Standard_Boolean& Empty,
-                        Standard_Boolean& Same,
-                        Standard_Boolean& Multpoint,
-                        IntPatch_SequenceOfLine& slin,
-                        IntPatch_SequenceOfPoint& spnt)
-
+Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
+                                          const IntSurf_Quadric& Quad2,
+                                          const IntAna_QuadQuadGeo& theInter,
+                                          const Standard_Real Tol,
+                                          const Standard_Boolean isTheReverse,
+                                          Standard_Boolean& Empty,
+                                          Standard_Boolean& Same,
+                                          Standard_Boolean& Multpoint,
+                                          IntPatch_SequenceOfLine& slin,
+                                          IntPatch_SequenceOfPoint& spnt)
 {
   IntPatch_Point ptsol;
-
+  
   Standard_Integer i;
 
   IntSurf_TypeTrans trans1,trans2;
@@ -453,39 +860,43 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
   gp_Cylinder Cy1(Quad1.Cylinder());
   gp_Cylinder Cy2(Quad2.Cylinder());
 
-  IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
-
-  if (!inter.IsDone())
-  {
-    return Standard_False;
-  }
-
-  typint = inter.TypeInter();
-  Standard_Integer NbSol = inter.NbSolutions();
+  typint = theInter.TypeInter();
+  Standard_Integer NbSol = theInter.NbSolutions();
   Empty = Standard_False;
   Same  = Standard_False;
 
   switch (typint)
-  {
+      {
   case IntAna_Empty:
     {
       Empty = Standard_True;
-    }
+      }
     break;
 
   case IntAna_Same:
-    {
+      {
       Same  = Standard_True;
-    }
+      }
     break;
-    
+
   case IntAna_Point:
     {
-      gp_Pnt psol(inter.Point(1));
-      Standard_Real U1,V1,U2,V2;
-      Quad1.Parameters(psol,U1,V1);
-      Quad2.Parameters(psol,U2,V2);
+      gp_Pnt psol(theInter.Point(1));
       ptsol.SetValue(psol,Tol,Standard_True);
+
+      Standard_Real U1,V1,U2,V2;
+
+      if(isTheReverse)
+      {
+        Quad2.Parameters(psol, U1, V1);
+        Quad1.Parameters(psol, U2, V2);
+      }
+      else
+      {
+        Quad1.Parameters(psol, U1, V1);
+        Quad2.Parameters(psol, U2, V2);
+      }
+
       ptsol.SetParameters(U1,V1,U2,V2);
       spnt.Append(ptsol);
     }
@@ -496,10 +907,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       gp_Pnt ptref;
       if (NbSol == 1)
       { // Cylinders are tangent to each other by line
-        linsol = inter.Line(1);
+        linsol = theInter.Line(1);
         ptref = linsol.Location();
+        
+        //Radius-vectors
         gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
         gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
+
+        //outer normal lines
         gp_Vec norm1(Quad1.Normale(ptref));
         gp_Vec norm2(Quad2.Normale(ptref));
         IntSurf_Situation situcyl1;
@@ -507,6 +922,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
 
         if (crb1.Dot(crb2) < 0.)
         { // centre de courbures "opposes"
+            //ATTENTION!!! 
+            //        Normal and Radius-vector of the 1st(!) cylinder
+            //        is used for judging what the situation of the 2nd(!)
+            //        cylinder is.
+
           if (norm1.Dot(crb1) > 0.)
           {
             situcyl2 = IntSurf_Inside;
@@ -526,7 +946,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
           }
         }
         else
-        {
+          {
           if (Cy1.Radius() < Cy2.Radius())
           {
             if (norm1.Dot(crb1) > 0.)
@@ -576,9 +996,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       {
         for (i=1; i <= NbSol; i++)
         {
-          linsol = inter.Line(i);
+          linsol = theInter.Line(i);
           ptref = linsol.Location();
           gp_Vec lsd = linsol.Direction();
+
+          //Theoretically, qwe = +/- 1.0.
           Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
           if (qwe >0.00000001)
           {
@@ -607,31 +1029,52 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       gp_Vec Tgt;
       gp_Pnt ptref;
       IntPatch_Point pmult1, pmult2;
-      
-      elipsol = inter.Ellipse(1);
-      
+
+      elipsol = theInter.Ellipse(1);
+
       gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
       gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
-      
+
       Multpoint = Standard_True;
       pmult1.SetValue(pttang1,Tol,Standard_True);
       pmult2.SetValue(pttang2,Tol,Standard_True);
       pmult1.SetMultiple(Standard_True);
       pmult2.SetMultiple(Standard_True);
-      
+
       Standard_Real oU1,oV1,oU2,oV2;
-      Quad1.Parameters(pttang1,oU1,oV1); 
-      Quad2.Parameters(pttang1,oU2,oV2);
+
+      if(isTheReverse)
+      {
+        Quad2.Parameters(pttang1,oU1,oV1);
+        Quad1.Parameters(pttang1,oU2,oV2);
+      }
+      else
+      {
+        Quad1.Parameters(pttang1,oU1,oV1);
+        Quad2.Parameters(pttang1,oU2,oV2);
+      }
+
       pmult1.SetParameters(oU1,oV1,oU2,oV2);
 
-      Quad1.Parameters(pttang2,oU1,oV1); 
-      Quad2.Parameters(pttang2,oU2,oV2);
+      if(isTheReverse)
+      {
+        Quad2.Parameters(pttang2,oU1,oV1); 
+        Quad1.Parameters(pttang2,oU2,oV2);        
+      }
+      else
+      {
+        Quad1.Parameters(pttang2,oU1,oV1); 
+        Quad2.Parameters(pttang2,oU2,oV2);
+      }
+
       pmult2.SetParameters(oU1,oV1,oU2,oV2);
-      
-      // on traite la premiere ellipse
 
+      // on traite la premiere ellipse
+        
       //-- Calcul de la Transition de la ligne 
       ElCLib::D1(0.,elipsol,ptref,Tgt);
+
+      //Theoretically, qwe = +/- |Tgt|.
       Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
       if (qwe>0.00000001)
       {
@@ -644,7 +1087,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         trans2 = IntSurf_Out;
       }
       else
-      { 
+      {
         trans1=trans2=IntSurf_Undecided; 
       }
 
@@ -660,8 +1103,18 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         aIP.SetValue(aP,Tol,Standard_False);
         aIP.SetMultiple(Standard_False);
         //
-        Quad1.Parameters(aP, aU1, aV1); 
-        Quad2.Parameters(aP, aU2, aV2);
+
+        if(isTheReverse)
+        {
+          Quad2.Parameters(aP, aU1, aV1); 
+          Quad1.Parameters(aP, aU2, aV2);          
+        }
+        else
+        {
+          Quad1.Parameters(aP, aU1, aV1); 
+          Quad2.Parameters(aP, aU2, aV2);
+        }
+
         aIP.SetParameters(aU1, aV1, aU2, aV2);
         //
         aIP.SetParameter(0.);
@@ -685,7 +1138,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       //-- Transitions calculee au point 0    OK
       //
       // on traite la deuxieme ellipse
-      elipsol = inter.Ellipse(2);
+      elipsol = theInter.Ellipse(2);
 
       Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
       Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
@@ -704,6 +1157,8 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       
       //-- Calcul des transitions de ligne pour la premiere ligne 
       ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);      
+
+      //Theoretically, qwe = +/- |Tgt|.
       qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
       if(qwe> 0.00000001)
       {
@@ -731,443 +1186,116 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         aIP.SetValue(aP,Tol,Standard_False);
         aIP.SetMultiple(Standard_False);
         //
-        Quad1.Parameters(aP, aU1, aV1); 
-        Quad2.Parameters(aP, aU2, aV2);
-        aIP.SetParameters(aU1, aV1, aU2, aV2);
-        //
-        aIP.SetParameter(0.);
-        glig->AddVertex(aIP);
-        glig->SetFirstPoint(1);
-        //
-        aIP.SetParameter(2.*M_PI);
-        glig->AddVertex(aIP);
-        glig->SetLastPoint(2);
-      }
-      //
-      glig->AddVertex(pmult1);
-      glig->AddVertex(pmult2);
-      //
-      slin.Append(glig);
-    }
-    break;
 
-  case IntAna_NoGeometricSolution:
-    {
-      Standard_Boolean bReverse;
-      Standard_Real U1,V1,U2,V2;
-      gp_Pnt psol;
-      //
-      bReverse=IsToReverse(Cy1, Cy2, Tol);
-      if (bReverse)
-      {
-        Cy2=Quad1.Cylinder();
-        Cy1=Quad2.Cylinder();
-      }
-      //
-      IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol);
-      if (!anaint.IsDone())
-      {
-        return Standard_False;
-      }
-      
-      if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0)
-      {
-        Empty = Standard_True;
-      }
-      else
-      {
-        NbSol = anaint.NbPnt();
-        for (i = 1; i <= NbSol; i++)
+        if(isTheReverse)
         {
-          psol = anaint.Point(i);
-          Quad1.Parameters(psol,U1,V1);
-          Quad2.Parameters(psol,U2,V2);
-          ptsol.SetValue(psol,Tol,Standard_True);
-          ptsol.SetParameters(U1,V1,U2,V2);
-          spnt.Append(ptsol);
+          Quad2.Parameters(aP, aU1, aV1); 
+          Quad1.Parameters(aP, aU2, aV2);          
         }
-
-        gp_Pnt ptvalid, ptf, ptl;
-        gp_Vec tgvalid;
-
-        Standard_Real first,last,para;
-        IntAna_Curve curvsol;
-        Standard_Boolean tgfound;
-        Standard_Integer kount;
-
-        NbSol = anaint.NbCurve();
-        for (i = 1; i <= NbSol; i++)
+        else
         {
-          curvsol = anaint.Curve(i);
-          curvsol.Domain(first,last);
-          ptf = curvsol.Value(first);
-          ptl = curvsol.Value(last);
-
-          para = last;
-          kount = 1;
-          tgfound = Standard_False;
-
-          while (!tgfound)
-          {
-            para = (1.123*first + para)/2.123;
-            tgfound = curvsol.D1u(para,ptvalid,tgvalid);
-            if (!tgfound)
-            {
-              kount ++;
-              tgfound = kount > 5;
-            }
-          }
-
-          Handle(IntPatch_ALine) alig;
-          if (kount <= 5)
-          {
-            Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid),
-                                                  Quad1.Normale(ptvalid));
-            if(qwe>0.00000001)
-            { 
-              trans1 = IntSurf_Out;
-              trans2 = IntSurf_In;
-            }
-            else if(qwe<-0.00000001)
-            {
-              trans1 = IntSurf_In;
-              trans2 = IntSurf_Out;
-            }
-            else
-            { 
-              trans1=trans2=IntSurf_Undecided; 
-            }
-            alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
-          }
-          else
-          {
-            alig = new IntPatch_ALine(curvsol,Standard_False);
-            //-- cout << "Transition indeterminee" << endl;
-          }
-
-          Standard_Boolean TempFalse1 = Standard_False;
-          Standard_Boolean TempFalse2 = Standard_False;
-
-          ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first,
-                        TempFalse2,ptl,last,Multpoint,Tol);
-          slin.Append(alig);
+          Quad1.Parameters(aP, aU1, aV1); 
+          Quad2.Parameters(aP, aU2, aV2);
         }
-      }
-    }
-    break;
-
-  default:
-    return Standard_False;
-  }
-
-  return Standard_True;
-}
-
-//=======================================================================
-//function : ShortCosForm
-//purpose  : Represents theCosFactor*cosA+theSinFactor*sinA as
-//            theCoeff*cos(A-theAngle) if it is possibly (all angles are
-//            in radians).
-//=======================================================================
-static void ShortCosForm( const Standard_Real theCosFactor,
-                          const Standard_Real theSinFactor,
-                          Standard_Real& theCoeff,
-                          Standard_Real& theAngle)
-{
-  theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
-  theAngle = 0.0;
-  if(IsEqual(theCoeff, 0.0))
-  {
-    theAngle = 0.0;
-    return;
-  }
-
-  theAngle = acos(Abs(theCosFactor/theCoeff));
 
-  if(theSinFactor > 0.0)
-  {
-    if(IsEqual(theCosFactor, 0.0))
-    {
-      theAngle = M_PI/2.0;
-    }
-    else if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI-theAngle;
-    }
-  }
-  else if(IsEqual(theSinFactor, 0.0))
-  {
-    if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI;
-    }
-  }
-  if(theSinFactor < 0.0)
-  {
-    if(theCosFactor > 0.0)
-    {
-      theAngle = 2.0*M_PI-theAngle;
-    }
-    else if(IsEqual(theCosFactor, 0.0))
-    {
-      theAngle = 3.0*M_PI/2.0;
-    }
-    else if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI+theAngle;
-    }
-  }
-}
-
-enum SearchBoundType
-{
-  SearchNONE = 0,
-  SearchV1 = 1,
-  SearchV2 = 2
-};
-
-//Stores equations coefficients
-struct stCoeffsValue
-{
-  stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
-
-  math_Vector mVecA1;
-  math_Vector mVecA2;
-  math_Vector mVecB1;
-  math_Vector mVecB2;
-  math_Vector mVecC1;
-  math_Vector mVecC2;
-  math_Vector mVecD;
-
-  Standard_Real mK21; //sinU2
-  Standard_Real mK11; //sinU1
-  Standard_Real mL21; //cosU2
-  Standard_Real mL11;  //cosU1
-  Standard_Real mM1;  //Free member
-
-  Standard_Real mK22; //sinU2
-  Standard_Real mK12; //sinU1
-  Standard_Real mL22; //cosU2
-  Standard_Real mL12; //cosU1
-  Standard_Real mM2; //Free member
-  
-  Standard_Real mK1;
-  Standard_Real mL1;
-  Standard_Real mK2;
-  Standard_Real mL2;
-
-  Standard_Real mFIV1;
-  Standard_Real mPSIV1;
-  Standard_Real mFIV2;
-  Standard_Real mPSIV2;
-
-  Standard_Real mB;
-  Standard_Real mC;
-  Standard_Real mFI1;
-  Standard_Real mFI2;
-};
-
-stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1,
-                              const gp_Cylinder& theCyl2):
-  mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
-  mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
-  mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
-  mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
-  mVecC1(theCyl1.Axis().Direction().XYZ()),
-  mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
-  mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
-{
-  enum CoupleOfEquation
-  {
-    COENONE = 0,
-    COE12 = 1,
-    COE23 = 2,
-    COE13 = 3
-  }aFoundCouple = COENONE;
-
-
-  Standard_Real aDetV1V2 = 0.0;
-  
-  const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
-  const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
-  const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
-  const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
-  const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
-  const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
-
-  if(anAbsD1 >= anAbsD2)
-  {
-    if(anAbsD3 > anAbsD1)
-    {
-      aFoundCouple = COE13;
-      aDetV1V2 = aDelta3;
-    }
-    else
-    {
-      aFoundCouple = COE12;
-      aDetV1V2 = aDelta1;
-    }
-  }
-  else
-  {
-    if(anAbsD3 > anAbsD2)
-    {
-      aFoundCouple = COE13;
-      aDetV1V2 = aDelta3;
-    }
-    else
-    {
-      aFoundCouple = COE23;
-      aDetV1V2 = aDelta2;
-    }
-  }
-
-  // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
-  // cross-product between directions (i.e. sine of angle).
-  // If sine is too small then sine is (approx.) equal to angle itself.
-  // Therefore, in this case we should compare sine with angular tolerance. 
-  // This constant is used for check if axes are parallel (see constructor
-  // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
-  if(Abs(aDetV1V2) < Precision::Angular())
-  {
-    Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
-  }
-
-  switch(aFoundCouple)
-  {
-  case COE12:
-    break;
-  case COE23:
-    {
-      math_Vector aVTemp(mVecA1);
-      mVecA1(1) = aVTemp(2);
-      mVecA1(2) = aVTemp(3);
-      mVecA1(3) = aVTemp(1);
-
-      aVTemp = mVecA2;
-      mVecA2(1) = aVTemp(2);
-      mVecA2(2) = aVTemp(3);
-      mVecA2(3) = aVTemp(1);
-
-      aVTemp = mVecB1;
-      mVecB1(1) = aVTemp(2);
-      mVecB1(2) = aVTemp(3);
-      mVecB1(3) = aVTemp(1);
-      
-      aVTemp = mVecB2;
-      mVecB2(1) = aVTemp(2);
-      mVecB2(2) = aVTemp(3);
-      mVecB2(3) = aVTemp(1);
-
-      aVTemp = mVecC1;
-      mVecC1(1) = aVTemp(2);
-      mVecC1(2) = aVTemp(3);
-      mVecC1(3) = aVTemp(1);
-
-      aVTemp = mVecC2;
-      mVecC2(1) = aVTemp(2);
-      mVecC2(2) = aVTemp(3);
-      mVecC2(3) = aVTemp(1);
-
-      aVTemp = mVecD;
-      mVecD(1) = aVTemp(2);
-      mVecD(2) = aVTemp(3);
-      mVecD(3) = aVTemp(1);
-
-      break;
-    }
-  case COE13:
-    {
-      math_Vector aVTemp = mVecA1;
-      mVecA1(2) = aVTemp(3);
-      mVecA1(3) = aVTemp(2);
-
-      aVTemp = mVecA2;
-      mVecA2(2) = aVTemp(3);
-      mVecA2(3) = aVTemp(2);
-
-      aVTemp = mVecB1;
-      mVecB1(2) = aVTemp(3);
-      mVecB1(3) = aVTemp(2);
-
-      aVTemp = mVecB2;
-      mVecB2(2) = aVTemp(3);
-      mVecB2(3) = aVTemp(2);
-
-      aVTemp = mVecC1;
-      mVecC1(2) = aVTemp(3);
-      mVecC1(3) = aVTemp(2);
-
-      aVTemp = mVecC2;
-      mVecC2(2) = aVTemp(3);
-      mVecC2(3) = aVTemp(2);
-
-      aVTemp = mVecD;
-      mVecD(2) = aVTemp(3);
-      mVecD(3) = aVTemp(2);
-
-      break;
+        aIP.SetParameters(aU1, aV1, aU2, aV2);
+        //
+        aIP.SetParameter(0.);
+        glig->AddVertex(aIP);
+        glig->SetFirstPoint(1);
+        //
+        aIP.SetParameter(2.*M_PI);
+        glig->AddVertex(aIP);
+        glig->SetLastPoint(2);
+      }
+      //
+      glig->AddVertex(pmult1);
+      glig->AddVertex(pmult2);
+      //
+      slin.Append(glig);
     }
-  default:
     break;
-  }
-
-  //------- For V1 (begin)
-  //sinU2
-  mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
-  //sinU1
-  mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
-  //cosU2
-  mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
-  //cosU1
-  mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
-  //Free member
-  mM1 =  (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
-  //------- For V1 (end)
-
-  //------- For V2 (begin)
-  //sinU2
-  mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
-  //sinU1
-  mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
-  //cosU2
-  mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
-  //cosU1
-  mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
-  //Free member
-  mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
-  //------- For V1 (end)
 
-  ShortCosForm(mL11, mK11, mK1, mFIV1);
-  ShortCosForm(mL21, mK21, mL1, mPSIV1);
-  ShortCosForm(mL12, mK12, mK2, mFIV2);
-  ShortCosForm(mL22, mK22, mL2, mPSIV2);
+  case IntAna_Parabola:
+  case IntAna_Hyperbola:
+    Standard_Failure::Raise("IntCyCy(): Wrong intersection type!");
 
-  const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
-                      aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
-                      aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
-                      aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
+  case IntAna_Circle:
+    // Circle is useful when we will work with trimmed surfaces
+    // (two cylinders can be tangent by their basises, e.g. circle)
+  case IntAna_NoGeometricSolution:
+  default:
+    return Standard_False;
+  }
 
-  mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2;  //Free
+  return Standard_True;
+}
 
-  Standard_Real aA = 0.0;
+//=======================================================================
+//function : ShortCosForm
+//purpose  : Represents theCosFactor*cosA+theSinFactor*sinA as
+//            theCoeff*cos(A-theAngle) if it is possibly (all angles are
+//            in radians).
+//=======================================================================
+static void ShortCosForm( const Standard_Real theCosFactor,
+                          const Standard_Real theSinFactor,
+                          Standard_Real& theCoeff,
+                          Standard_Real& theAngle)
+{
+  theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor);
+  theAngle = 0.0;
+  if(IsEqual(theCoeff, 0.0))
+  {
+    theAngle = 0.0;
+    return;
+  }
 
-  ShortCosForm(aB2,aB1,mB,mFI1);
-  ShortCosForm(aA2,aA1,aA,mFI2);
+  theAngle = acos(Abs(theCosFactor/theCoeff));
 
-  mB /= aA;
-  mC /= aA;
+  if(theSinFactor > 0.0)
+  {
+    if(IsEqual(theCosFactor, 0.0))
+    {
+      theAngle = M_PI/2.0;
+    }
+    else if(theCosFactor < 0.0)
+    {
+      theAngle = M_PI-theAngle;
+    }
+  }
+  else if(IsEqual(theSinFactor, 0.0))
+  {
+    if(theCosFactor < 0.0)
+    {
+      theAngle = M_PI;
+    }
+  }
+  if(theSinFactor < 0.0)
+  {
+    if(theCosFactor > 0.0)
+    {
+      theAngle = 2.0*M_PI-theAngle;
+    }
+    else if(IsEqual(theCosFactor, 0.0))
+    {
+      theAngle = 3.0*M_PI/2.0;
+    }
+    else if(theCosFactor < 0.0)
+    {
+      theAngle = M_PI+theAngle;
+    }
+  }
 }
 
 //=======================================================================
 //function : CylCylMonotonicity
 //purpose  : Determines, if U2(U1) function is increasing.
 //=======================================================================
-static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par,
-                                            const Standard_Integer theWLIndex,
-                                            const stCoeffsValue& theCoeffs,
-                                            const Standard_Real thePeriod,
-                                            Standard_Boolean& theIsIncreasing)
+Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
+                                                        const Standard_Integer theWLIndex,
+                                                        const stCoeffsValue& theCoeffs,
+                                                        const Standard_Real thePeriod,
+                                                        Standard_Boolean& theIsIncreasing)
 {
   // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
   //Accordingly,
@@ -1235,11 +1363,11 @@ static Standard_Boolean CylCylMonotonicity( const Standard_Real theU1par,
 //            esimates the tolerance of U2-computing (estimation result is
 //            assigned to *theDelta value).
 //=======================================================================
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
                                                 const Standard_Integer theWLIndex,
                                                 const stCoeffsValue& theCoeffs,
                                                 Standard_Real& theU2,
-                                                Standard_Real* const theDelta = 0)
+                                                Standard_Real* const theDelta)
 {
   //This formula is got from some experience and can be changed.
   const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
@@ -1308,7 +1436,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
   return Standard_True;
 }
 
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
                                                 const Standard_Real theU2,
                                                 const stCoeffsValue& theCoeffs,
                                                 Standard_Real& theV1,
@@ -1328,7 +1456,7 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
 }
 
 
-static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
                                                 const Standard_Integer theWLIndex,
                                                 const stCoeffsValue& theCoeffs,
                                                 Standard_Real& theU2,
@@ -1348,13 +1476,13 @@ static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
 //function : SearchOnVBounds
 //purpose  : 
 //=======================================================================
-static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
-                                        const stCoeffsValue& theCoeffs,
+Standard_Boolean WorkWithBoundaries::
+                        SearchOnVBounds(const SearchBoundType theSBType,
                                         const Standard_Real theVzad,
                                         const Standard_Real theVInit,
                                         const Standard_Real theInitU2,
                                         const Standard_Real theInitMainVar,
-                                        Standard_Real& theMainVariableValue)
+                                        Standard_Real& theMainVariableValue) const
 {
   const Standard_Integer aNbDim = 3;
   const Standard_Real aMaxError = 4.0*M_PI; // two periods
@@ -1381,40 +1509,40 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
                         aSinU2 = sin(aU2Prev),
                         aCosU2 = cos(aU2Prev);
 
-    math_Vector aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev +
-                                              theCoeffs.mVecB2) * aSinU2 -
-                              (theCoeffs.mVecB2 * aU2Prev -
-                                              theCoeffs.mVecA2) * aCosU2 +
-                              (theCoeffs.mVecA1 * aMainVarPrev +
-                                              theCoeffs.mVecB1) * aSinU1 -
-                              (theCoeffs.mVecB1 * aMainVarPrev -
-                                              theCoeffs.mVecA1) * aCosU1 +
-                                                            theCoeffs.mVecD;
+    math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
+                                              myCoeffs.mVecB2) * aSinU2 -
+                              (myCoeffs.mVecB2 * aU2Prev -
+                                              myCoeffs.mVecA2) * aCosU2 +
+                              (myCoeffs.mVecA1 * aMainVarPrev +
+                                              myCoeffs.mVecB1) * aSinU1 -
+                              (myCoeffs.mVecB1 * aMainVarPrev -
+                                              myCoeffs.mVecA1) * aCosU1 +
+                                                            myCoeffs.mVecD;
 
     math_Vector aMSum(1, 3);
 
     switch(theSBType)
     {
     case SearchV1:
-      aMatr.SetCol(3, theCoeffs.mVecC2);
-      aMSum = theCoeffs.mVecC1 * theVzad;
+      aMatr.SetCol(3, myCoeffs.mVecC2);
+      aMSum = myCoeffs.mVecC1 * theVzad;
       aVecFreeMem -= aMSum;
-      aMSum += theCoeffs.mVecC2*anOtherVar;
+      aMSum += myCoeffs.mVecC2*anOtherVar;
       break;
 
     case SearchV2:
-      aMatr.SetCol(3, theCoeffs.mVecC1);
-      aMSum = theCoeffs.mVecC2 * theVzad;
+      aMatr.SetCol(3, myCoeffs.mVecC1);
+      aMSum = myCoeffs.mVecC2 * theVzad;
       aVecFreeMem -= aMSum;
-      aMSum += theCoeffs.mVecC1*anOtherVar;
+      aMSum += myCoeffs.mVecC1*anOtherVar;
       break;
 
     default:
       return Standard_False;
     }
 
-    aMatr.SetCol(1, theCoeffs.mVecA1 * aSinU1 - theCoeffs.mVecB1 * aCosU1);
-    aMatr.SetCol(2, theCoeffs.mVecA2 * aSinU2 - theCoeffs.mVecB2 * aCosU2);
+    aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
+    aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
 
     Standard_Real aDetMainSyst = aMatr.Determinant();
 
@@ -1460,11 +1588,11 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
                           aCosU1Last = cos(aMainVarPrev),
                           aSinU2Last = sin(aU2Prev),
                           aCosU2Last = cos(aU2Prev);
-      aMSum -= (theCoeffs.mVecA1*aCosU1Last + 
-                theCoeffs.mVecB1*aSinU1Last +
-                theCoeffs.mVecA2*aCosU2Last +
-                theCoeffs.mVecB2*aSinU2Last +
-                theCoeffs.mVecD);
+      aMSum -= (myCoeffs.mVecA1*aCosU1Last + 
+                myCoeffs.mVecB1*aSinU1Last +
+                myCoeffs.mVecA2*aCosU2Last +
+                myCoeffs.mVecB2*aSinU2Last +
+                myCoeffs.mVecD);
       const Standard_Real aSQNorm = aMSum.Norm2();
       return (aSQNorm < aTol2);
     }
@@ -1531,27 +1659,10 @@ Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
     return Standard_True;
   }
 
-  if(IsEqual(thePeriod, 0.0))
-  {//it cannot be inscribed
-    return Standard_False;
-  }
-
-  //Make theUGiven nearer to theUfTarget (in order to avoid
-  //excess iterations)
-  theUGiven += thePeriod*Floor((theUfTarget-theUGiven)/thePeriod);
-  Standard_Real aDU = theUGiven - theUfTarget;
-
-  if(aDU > 0.0)
-    aDU = -thePeriod;
-  else
-    aDU = +thePeriod;
+  const Standard_Real aUf = theUfTarget - theTol2D;
+  const Standard_Real aUl = aUf + thePeriod;
 
-  while(((theUGiven - theUfTarget)*aDU < 0.0) && 
-        !((theUfTarget - theUGiven <= theTol2D) &&
-        (theUGiven - theUlTarget <= theTol2D)))
-  {
-    theUGiven += aDU;
-  }
+  theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
   
   return ((theUfTarget - theUGiven <= theTol2D) &&
           (theUGiven - theUlTarget <= theTol2D));
@@ -1625,7 +1736,7 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
     if((anA-anB) < theTol)
     {
       if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l)) 
-        anA = (anA + anB)/2.0;
+      anA = (anA + anB)/2.0;
       else
         anA = anB;
 
@@ -1645,7 +1756,7 @@ static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
 //=======================================================================
 static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
                                         const IntSurf_Quadric& theQuad2,
-                                        const stCoeffsValue& theCoeffs,
+                                        const ComputationMethods::stCoeffsValue& theCoeffs,
                                         const Standard_Boolean isTheReverse,
                                         const Standard_Boolean isThePrecise,
                                         const gp_Pnt2d& thePntOnSurf1,
@@ -1713,8 +1824,21 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
       aPlast.ParametersOnS1(aUl, aVl);
 
     if(!theFlBefore && (aU1par <= aUl))
-    {//Parameter value must be increased if theFlBefore == FALSE.
-      return Standard_False;
+    {
+      //Parameter value must be increased if theFlBefore == FALSE.
+      
+      aU1par += thePeriodOfSurf1;
+
+      //The condition is as same as in
+      //InscribePoint(...) function
+      if((theUfSurf1 - aU1par > theTol2D) ||
+         (aU1par - theUlSurf1 > theTol2D))
+      {
+        //New aU1par is out of target interval.
+        //Go back to old value.
+
+        return Standard_False;
+      }
     }
 
     //theTol2D is minimal step along parameter changed. 
@@ -1760,8 +1884,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
     {
       SeekAdditionalPoints( theQuad1, theQuad2, theLine, 
                             theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
-                            aNbPnts-1, theUfSurf2, theUlSurf2,
-                            theTol2D, thePeriodOfSurf1, isTheReverse);
+                            aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
     }
   }
 
@@ -1772,15 +1895,7 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
 //function : AddBoundaryPoint
 //purpose  : 
 //=======================================================================
-static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
-                                          const IntSurf_Quadric& theQuad2,
-                                          const Handle(IntPatch_WLine)& theWL,
-                                          const stCoeffsValue& theCoeffs,
-                                          const Bnd_Box2d& theUVSurf1,
-                                          const Bnd_Box2d& theUVSurf2,
-                                          const Standard_Real theTol3D,
-                                          const Standard_Real theTol2D,
-                                          const Standard_Real thePeriod,
+void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
                                           const Standard_Real theU1,
                                           const Standard_Real theU2,
                                           const Standard_Real theV1,
@@ -1788,10 +1903,9 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
                                           const Standard_Real theV2,
                                           const Standard_Real theV2Prev,
                                           const Standard_Integer theWLIndex,
-                                          const Standard_Boolean isTheReverse,
                                           const Standard_Boolean theFlForce,
                                           Standard_Boolean& isTheFound1,
-                                          Standard_Boolean& isTheFound2)
+                                          Standard_Boolean& isTheFound2) const
 {
   Standard_Real aUSurf1f = 0.0, //const
                 aUSurf1l = 0.0,
@@ -1802,176 +1916,97 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
                 aVSurf2f = 0.0,
                 aVSurf2l = 0.0;
 
-  theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
-  theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
-
-  SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
-  Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
+  myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
+  myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
+  
+  const Standard_Integer aSize = 4;
+  const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
 
-  Standard_Real anUpar1 = theU1, anUpar2 = theU1;
-  Standard_Real aVf = theV1, aVl = theV1Prev;
+  StPInfo aUVPoint[aSize];
 
-  if( (Abs(aVf-aVSurf1f) <= theTol2D) ||
-      ((aVf-aVSurf1f)*(aVl-aVSurf1f) <= 0.0))
-  {
-    aTS1 = SearchV1;
-    aV1zad = aVSurf1f;
-    isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theV2, theU2, theU1, anUpar1);
-  }
-  else if((Abs(aVf-aVSurf1l) <= theTol2D) ||
-          ((aVf-aVSurf1l)*(aVl-aVSurf1l) <= 0.0))
+  for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
   {
-    aTS1 = SearchV1;
-    aV1zad = aVSurf1l;
-    isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theV2, theU2, theU1, anUpar1);
-  }
+    const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
+                        aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
 
-  aVf = theV2;
-  aVl = theV2Prev;
+    const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
 
-  if((Abs(aVf-aVSurf2f) <= theTol2D) || 
-      ((aVf-aVSurf2f)*(aVl-aVSurf2f) <= 0.0))
-  {
-    aTS2 = SearchV2;
-    aV2zad = aVSurf2f;
-    isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theV1, theU2, theU1, anUpar2);
-  }
-  else if((Abs(aVf-aVSurf2l) <= theTol2D) ||
-          ((aVf-aVSurf2l)*(aVl-aVSurf2l) <= 0.0))
-  {
-    aTS2 = SearchV2;
-    aV2zad = aVSurf2l;
-    isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theV1, theU2, theU1, anUpar2);
-  }
+    for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
+    {
 
-  if(!isTheFound1 && !isTheFound2)
-    return Standard_True;
+      const Standard_Integer anIndex = anIDSurf+anIDBound;
 
-  //We increase U1 parameter. Therefore, added point must have U1 parameter less than
-  //or equal to current (conditions "(anUpar1 < theU1)" and "(anUpar2 < theU1)").
+      aUVPoint[anIndex].mySurfID = anIDSurf;
 
-  if(anUpar1 < anUpar2)
-  {
-    if(isTheFound1 && (anUpar1 < theU1))
-    {
-      Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
-      if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
-      {
-        isTheFound1 = Standard_False;
-      }
-      
-      if(aTS1 == SearchV1)
-        aV1 = aV1zad;
-
-      if(aTS1 == SearchV2)
-        aV2 = aV2zad;
-
-      if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
-                                        gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
-                                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
-                                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
-                                        theWL->Curve(), theWLIndex, theTol3D,
-                                        theTol2D, theFlForce))
+      if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
+          ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
       {
-        isTheFound1 = Standard_False;
+        continue;
       }
-    }
-    else
-    {
-      isTheFound1 = Standard_False;
-    }
 
-    if(isTheFound2 && (anUpar2 < theU1))
-    {
-      Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
-      if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
+      const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
+                                                    (anIDSurf == 0) ? theV2 : theV1,
+                                                    theU2, theU1,
+                                                    aUVPoint[anIndex].myU1);
+
+      if(!aRes || aUVPoint[anIndex].myU1 >= theU1)
       {
-        isTheFound2 = Standard_False;
+        aUVPoint[anIndex].myU1 = RealLast();
+        continue;
       }
-
-      if(aTS2 == SearchV1)
-        aV1 = aV1zad;
-
-      if(aTS2 == SearchV2)
-        aV2 = aV2zad;
-
-      if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
-                                        gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
-                                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
-                                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
-                                        theWL->Curve(), theWLIndex, theTol3D,
-                                        theTol2D, theFlForce))
+      else
       {
-        isTheFound2 = Standard_False;
+        Standard_Real &aU1 = aUVPoint[anIndex].myU1,
+                      &aU2 = aUVPoint[anIndex].myU2,
+                      &aV1 = aUVPoint[anIndex].myV1,
+                      &aV2 = aUVPoint[anIndex].myV2;
+        aU2 = theU2;
+        aV1 = theV1;
+        aV2 = theV2;
+
+        if(!ComputationMethods::
+                  CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
+        {
+          aU1 = RealLast();
+          continue;
+        }
+
+        if(aTS == SearchV1)
+          aV1 = anArrVzad[anIndex];
+        else //if(aTS[anIndex] == SearchV2)
+          aV2 = anArrVzad[anIndex];
       }
     }
-    else
-    {
-      isTheFound2 = Standard_False;
-    }
   }
-  else
-  {
-    if(isTheFound2 && (anUpar2 < theU1))
-    {
-      Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
-      if(!CylCylComputeParameters(anUpar2, theWLIndex, theCoeffs, aU2, aV1, aV2))
-      {
-        isTheFound2 = Standard_False;
-      }
 
-      if(aTS2 == SearchV1)
-        aV1 = aV1zad;
+  std::sort(aUVPoint, aUVPoint+aSize);
+    
+  isTheFound1 = isTheFound2 = Standard_False;
 
-      if(aTS2 == SearchV2)
-        aV2 = aV2zad;
+  for(Standard_Integer i = 0; i < aSize; i++)
+  {
+    if(aUVPoint[i].myU1 == RealLast())
+      break;
 
-      if(isTheFound2 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
-                                        gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
-                                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
-                                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
-                                        theWL->Curve(), theWLIndex, theTol3D,
-                                        theTol2D, theFlForce))
-      {
-        isTheFound2 = Standard_False;
-      }
-    }
-    else
+    if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
+                        gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
+                        gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
+                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
+                        theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
     {
-      isTheFound2 = Standard_False;
+      continue;
     }
 
-    if(isTheFound1 && (anUpar1 < theU1))
+    if(aUVPoint[i].mySurfID == 0)
     {
-      Standard_Real aU2 = theU2, aV1 = theV1, aV2 = theV2;
-      if(!CylCylComputeParameters(anUpar1, theWLIndex, theCoeffs, aU2, aV1, aV2))
-      {
-        isTheFound1 = Standard_False;
-      }
-
-      if(aTS1 == SearchV1)
-        aV1 = aV1zad;
-
-      if(aTS1 == SearchV2)
-        aV2 = aV2zad;
-
-      if(isTheFound1 && !AddPointIntoWL(theQuad1, theQuad2, theCoeffs, isTheReverse, Standard_False,
-                                        gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
-                                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
-                                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, thePeriod,
-                                        theWL->Curve(), theWLIndex, theTol3D,
-                                        theTol2D, theFlForce))
-      {
-        isTheFound1 = Standard_False;
-      }
+      isTheFound1 = Standard_True;
     }
     else
     {
-      isTheFound1 = Standard_False;
+      isTheFound2 = Standard_True;
     }
   }
-
-  return Standard_True;
 }
 
 //=======================================================================
@@ -1981,13 +2016,11 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const IntSurf_Quadric& theQuad2,
                                   const Handle(IntSurf_LineOn2S)& theLine,
-                                  const stCoeffsValue& theCoeffs,
+                                  const ComputationMethods::stCoeffsValue& theCoeffs,
                                   const Standard_Integer theWLIndex,
                                   const Standard_Integer theMinNbPoints,
                                   const Standard_Integer theStartPointOnLine,
                                   const Standard_Integer theEndPointOnLine,
-                                  const Standard_Real theU2f,
-                                  const Standard_Real theU2l,
                                   const Standard_Real theTol2D,
                                   const Standard_Real thePeriodOfSurf2,
                                   const Standard_Boolean isTheReverse)
@@ -2062,16 +2095,23 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
 
       U1prec = 0.5*(U1f+U1l);
       
-      if(!CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
+      if(!ComputationMethods::
+            CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
+      {
         continue;
+      }
 
-      InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
+      MinMax(U2f, U2l);
+      if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
+      {
+        continue;
+      }
 
       const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
       const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
 
-#ifdef OCCT_DEBUG
-      //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+      std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
 #endif
 
       IntSurf_PntOn2S anIP;
@@ -2102,27 +2142,24 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
 //purpose  : Computes true domain of future intersection curve (allows
 //            avoiding excess iterations)
 //=======================================================================
-//=======================================================================
-static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
-                                            const Standard_Real thePeriod,
-                                            Standard_Real theU1f[],
-                                            Standard_Real theU1l[])
+Standard_Boolean WorkWithBoundaries::BoundariesComputing(Standard_Real theU1f[],
+                                                         Standard_Real theU1l[]) const
 {
-  if(theCoeffs.mB > 0.0)
+  if(myCoeffs.mB > 0.0)
   {
-    if(theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
+    if(myCoeffs.mB + Abs(myCoeffs.mC) < -1.0)
     {//There is NOT intersection
       return Standard_False;
     }
-    else if(theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    else if(myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0)
     {//U=[0;2*PI]+aFI1
-      theU1f[0] = theCoeffs.mFI1;
-      theU1l[0] = thePeriod + theCoeffs.mFI1;
+      theU1f[0] = myCoeffs.mFI1;
+      theU1l[0] = myPeriod + myCoeffs.mFI1;
     }
-    else if((1 + theCoeffs.mC <= theCoeffs.mB) &&
-            (theCoeffs.mB <= 1 - theCoeffs.mC))
+    else if((1 + myCoeffs.mC <= myCoeffs.mB) &&
+            (myCoeffs.mB <= 1 - myCoeffs.mC))
     {
-      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
@@ -2130,15 +2167,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
 
       const Standard_Real aDAngle = acos(anArg);
       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
-      theU1f[0] = theCoeffs.mFI1;
-      theU1l[0] = aDAngle + theCoeffs.mFI1;
-      theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1;
-      theU1l[1] = thePeriod + theCoeffs.mFI1;
+      theU1f[0] = myCoeffs.mFI1;
+      theU1l[0] = aDAngle + myCoeffs.mFI1;
+      theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1;
+      theU1l[1] = myPeriod + myCoeffs.mFI1;
     }
-    else if((1 - theCoeffs.mC <= theCoeffs.mB) &&
-            (theCoeffs.mB <= 1 + theCoeffs.mC))
+    else if((1 - myCoeffs.mC <= myCoeffs.mB) &&
+            (myCoeffs.mB <= 1 + myCoeffs.mC))
     {
-      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
+      Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
@@ -2147,13 +2184,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
       const Standard_Real aDAngle = acos(anArg);
       //U=[aDAngle;2*PI-aDAngle]+aFI1
 
-      theU1f[0] = aDAngle + theCoeffs.mFI1;
-      theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1;
+      theU1f[0] = aDAngle + myCoeffs.mFI1;
+      theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1;
     }
-    else if(theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
+    else if(myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0)
     {
-      Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
-                    anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      Standard_Real anArg1 = (1 - myCoeffs.mC) / myCoeffs.mB,
+                    anArg2 = -(myCoeffs.mC + 1) / myCoeffs.mB;
       if(anArg1 > 1.0)
         anArg1 = 1.0;
       if(anArg1 < -1.0)
@@ -2168,31 +2205,31 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
       //(U=[aDAngle1;aDAngle2]+aFI1) ||
       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
 
-      theU1f[0] = aDAngle1 + theCoeffs.mFI1;
-      theU1l[0] = aDAngle2 + theCoeffs.mFI1;
-      theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1;
-      theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1;
+      theU1f[0] = aDAngle1 + myCoeffs.mFI1;
+      theU1l[0] = aDAngle2 + myCoeffs.mFI1;
+      theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1;
+      theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1;
     }
     else
     {
       return Standard_False;
     }
   }
-  else if(theCoeffs.mB < 0.0)
+  else if(myCoeffs.mB < 0.0)
   {
-    if(theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
+    if(myCoeffs.mB + Abs(myCoeffs.mC) > 1.0)
     {//There is NOT intersection
       return Standard_False;
     }
-    else if(-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    else if(-myCoeffs.mB + Abs(myCoeffs.mC) <= 1.0)
     {//U=[0;2*PI]+aFI1
-      theU1f[0] = theCoeffs.mFI1;
-      theU1l[0] = thePeriod + theCoeffs.mFI1;
+      theU1f[0] = myCoeffs.mFI1;
+      theU1l[0] = myPeriod + myCoeffs.mFI1;
     }
-    else if((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
-            ( theCoeffs.mB <= theCoeffs.mC - 1))
+    else if((-myCoeffs.mC - 1 <= myCoeffs.mB) &&
+            ( myCoeffs.mB <= myCoeffs.mC - 1))
     {
-      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
+      Standard_Real anArg = (1 - myCoeffs.mC) / myCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
@@ -2201,15 +2238,15 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
       const Standard_Real aDAngle = acos(anArg);
       //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
 
-      theU1f[0] = theCoeffs.mFI1;
-      theU1l[0] = aDAngle + theCoeffs.mFI1;
-      theU1f[1] = thePeriod - aDAngle + theCoeffs.mFI1;
-      theU1l[1] = thePeriod + theCoeffs.mFI1;
+      theU1f[0] = myCoeffs.mFI1;
+      theU1l[0] = aDAngle + myCoeffs.mFI1;
+      theU1f[1] = myPeriod - aDAngle + myCoeffs.mFI1;
+      theU1l[1] = myPeriod + myCoeffs.mFI1;
     }
-    else if((theCoeffs.mC - 1 <= theCoeffs.mB) &&
-            (theCoeffs.mB <= -theCoeffs.mB - 1))
+    else if((myCoeffs.mC - 1 <= myCoeffs.mB) &&
+            (myCoeffs.mB <= -myCoeffs.mB - 1))
     {
-      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
+      Standard_Real anArg = -(myCoeffs.mC + 1) / myCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
@@ -2218,13 +2255,13 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
       const Standard_Real aDAngle = acos(anArg);
       //U=[aDAngle;2*PI-aDAngle]+aFI1
 
-      theU1f[0] = aDAngle + theCoeffs.mFI1;
-      theU1l[0] = thePeriod - aDAngle + theCoeffs.mFI1;
+      theU1f[0] = aDAngle + myCoeffs.mFI1;
+      theU1l[0] = myPeriod - aDAngle + myCoeffs.mFI1;
     }
-    else if(-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
+    else if(-myCoeffs.mB - Abs(myCoeffs.mC) >= 1.0)
     {
-      Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
-                    anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
+      Standard_Real anArg1 = -(myCoeffs.mC + 1) / myCoeffs.mB,
+                    anArg2 = (1 - myCoeffs.mC) / myCoeffs.mB;
       if(anArg1 > 1.0)
         anArg1 = 1.0;
       if(anArg1 < -1.0)
@@ -2239,10 +2276,10 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
       //(U=[aDAngle1;aDAngle2]+aFI1) ||
       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
 
-      theU1f[0] = aDAngle1 + theCoeffs.mFI1;
-      theU1l[0] = aDAngle2 + theCoeffs.mFI1;
-      theU1f[1] = thePeriod - aDAngle2 + theCoeffs.mFI1;
-      theU1l[1] = thePeriod - aDAngle1 + theCoeffs.mFI1;
+      theU1f[0] = aDAngle1 + myCoeffs.mFI1;
+      theU1l[0] = aDAngle2 + myCoeffs.mFI1;
+      theU1f[1] = myPeriod - aDAngle2 + myCoeffs.mFI1;
+      theU1l[1] = myPeriod - aDAngle1 + myCoeffs.mFI1;
     }
     else
     {
@@ -2262,7 +2299,7 @@ static Standard_Boolean BoundariesComputing(const stCoeffsValue& theCoeffs,
 //purpose  : theNbCritPointsMax contains true number of critical points.
 //            It must be initialized correctly before function calling
 //=======================================================================
-static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
+static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
                                     const Standard_Real theUSurf1f,
                                     const Standard_Real theUSurf1l,
                                     const Standard_Real theUSurf2f,
@@ -2272,13 +2309,13 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
                                     Standard_Integer& theNbCritPointsMax,
                                     Standard_Real theU1crit[])
 {
-  //[0...1]   - in these points parameter U1 goes through
-  //            the seam-edge of the first cylinder.
-  //[2...3]   - First and last U1 parameter.
-  //[4...5]   - in these points parameter U2 goes through
-  //            the seam-edge of the second cylinder.
-  //[6...9]   - in these points an intersection line goes through
-  //            U-boundaries of the second surface.
+  //[0...1] - in these points parameter U1 goes through
+  //          the seam-edge of the first cylinder.
+  //[2...3] - First and last U1 parameter.
+  //[4...5] - in these points parameter U2 goes through
+  //          the seam-edge of the second cylinder.
+  //[6...9] - in these points an intersection line goes through
+  //          U-boundaries of the second surface.
   //[10...11] - Boundary of monotonicity interval of U2(U1) function
   //            (see CylCylMonotonicity() function)
 
@@ -2385,20 +2422,108 @@ static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
 }
 
 //=======================================================================
-//function : IntCyCyTrim
+//function : BoundaryEstimation
+//purpose  : Rough estimation of the parameter range.
+//=======================================================================
+void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
+                                            const gp_Cylinder& theCy2,
+                                            Bnd_Range& theOutBoxS1,
+                                            Bnd_Range& theOutBoxS2) const
+{
+  const gp_Dir &aD1 = theCy1.Axis().Direction(),
+               &aD2 = theCy2.Axis().Direction();
+  const Standard_Real aR1 = theCy1.Radius(),
+                      aR2 = theCy2.Radius();
+
+  //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
+  //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
+  //In fact, this parallelogram is a projection of the cylinders to the plane
+  //created by the intersected axes aD1 and aD2 (if the axes are skewed then
+  //one axis can be translated by parallel shifting till intersection).
+
+  const Standard_Real aCosA = aD1.Dot(aD2);
+  const Standard_Real aSqSinA = 1.0 - aCosA * aCosA;
+
+  //If sine is small then it can be compared with angle.
+  if (aSqSinA < Precision::Angular()*Precision::Angular())
+    return;
+
+  //Half of delta V. Delta V is a distance between 
+  //projections of two opposite parallelogram vertices
+  //(joined by the maximal diagonal) to the cylinder axis.
+  const Standard_Real aSinA = sqrt(aSqSinA);
+  const Standard_Real aHDV1 = (aR1 * aCosA + aR2)/aSinA,
+                      aHDV2 = (aR2 * aCosA + aR1)/aSinA;
+
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+  //The code in this block is created for test only.It is stupidly to create
+  //OCCT-test for the method, which will be changed possibly never.
+  std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
+#endif
+
+  //V-parameters of intersection point of the axes (in case of skewed axes,
+  //see comment above).
+  Standard_Real aV01 = 0.0, aV02 = 0.0;
+  ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);
+
+  theOutBoxS1.Add(aV01 - aHDV1);
+  theOutBoxS1.Add(aV01 + aHDV1);
+
+  theOutBoxS2.Add(aV02 - aHDV2);
+  theOutBoxS2.Add(aV02 + aHDV2);
+
+  theOutBoxS1.Enlarge(Precision::Confusion());
+  theOutBoxS2.Enlarge(Precision::Confusion());
+
+  Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
+  myUVSurf1.Get(aU1, aV1, aU2, aV2);
+  theOutBoxS1.Common(Bnd_Range(aV1, aV2));
+
+  myUVSurf2.Get(aU1, aV1, aU2, aV2);
+  theOutBoxS2.Common(Bnd_Range(aV1, aV2));
+}
+
+//=======================================================================
+//function : IntCyCy
 //purpose  : 
 //=======================================================================
-Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
-                              const IntSurf_Quadric& theQuad2,
-                              const Standard_Real theTol3D,
-                              const Standard_Real theTol2D,
-                              const Bnd_Box2d& theUVSurf1,
-                              const Bnd_Box2d& theUVSurf2,
-                              const Standard_Boolean isTheReverse,
-                              Standard_Boolean& isTheEmpty,
-                              IntPatch_SequenceOfLine& theSlin,
-                              IntPatch_SequenceOfPoint& theSPnt)
+Standard_Boolean IntCyCy( const IntSurf_Quadric& theQuad1,
+                          const IntSurf_Quadric& theQuad2,
+                          const Standard_Real theTol3D,
+                          const Standard_Real theTol2D,
+                          const Bnd_Box2d& theUVSurf1,
+                          const Bnd_Box2d& theUVSurf2,
+                          const Standard_Boolean isTheReverse,
+                          Standard_Boolean& isTheEmpty,
+                          Standard_Boolean& isTheSameSurface,
+                          Standard_Boolean& isTheMultiplePoint,
+                          IntPatch_SequenceOfLine& theSlin,
+                          IntPatch_SequenceOfPoint& theSPnt)
 {
+  isTheEmpty = Standard_True;
+  isTheSameSurface = Standard_False;
+  isTheMultiplePoint = Standard_False;
+  theSlin.Clear();
+  theSPnt.Clear();
+
+  const gp_Cylinder &aCyl1 = theQuad1.Cylinder(),
+                    &aCyl2 = theQuad2.Cylinder();
+
+  IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
+
+  if (!anInter.IsDone())
+  {
+    return Standard_False;
+  }
+
+  if(anInter.TypeInter() != IntAna_NoGeometricSolution)
+  {
+    return CyCyAnalyticalIntersect( theQuad1, theQuad2, anInter,
+                                    theTol3D, isTheReverse, isTheEmpty,
+                                    isTheSameSurface, isTheMultiplePoint,
+                                    theSlin, theSPnt);
+  }
+  
   Standard_Real aUSurf1f = 0.0, //const
                 aUSurf1l = 0.0,
                 aVSurf1f = 0.0,
@@ -2411,27 +2536,6 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
   theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
   theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
 
-  const gp_Cylinder&  aCyl1 = theQuad1.Cylinder(),
-                      aCyl2 = theQuad2.Cylinder();
-
-  IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
-
-  if (!anInter.IsDone())
-  {
-    return Standard_False;
-  }
-
-  IntAna_ResultType aTypInt = anInter.TypeInter();
-
-  if(aTypInt != IntAna_NoGeometricSolution)
-  { //It is not necessary (because result is an analytic curve) or
-    //it is impossible to make Walking-line.
-
-    return Standard_False;
-  }
-    
-  theSlin.Clear();
-  theSPnt.Clear();
   const Standard_Integer aNbMaxPoints = 2000;
   const Standard_Integer aNbMinPoints = 200;
   const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, 
@@ -2444,14 +2548,23 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
 
   const Standard_Integer aNbWLines = 2;
 
-  const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
+  const ComputationMethods::stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
+
+  const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs,
+                                      theUVSurf1, theUVSurf2, aNbWLines,
+                                      aPeriod, theTol3D, theTol2D, isTheReverse);
+
+  Bnd_Range aRangeS1, aRangeS2;
+  aBoundWork.BoundaryEstimation(aCyl1, aCyl2, aRangeS1, aRangeS2);
+  if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
+    return Standard_True;
 
   //Boundaries
   const Standard_Integer aNbOfBoundaries = 2;
   Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
   Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
 
-  if(!BoundariesComputing(anEquationCoeffs, aPeriod, aU1f, aU1l))
+  if(!aBoundWork.BoundariesComputing(aU1f, aU1l))
     return Standard_True;
 
   for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
@@ -2620,7 +2733,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
            //or on seam-edge strictly (if it is possible).
 
             Standard_Real aTol = theTol2D;
-            CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i], &aTol);
+            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
+                                                        aU2[i], &aTol);
             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
 
             aTol = Max(aTol, theTol2D);
@@ -2636,7 +2750,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
           }
           else
           {
-            CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
+            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
             InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
           }
           
@@ -2651,7 +2765,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
               //correct value.
 
               Standard_Boolean isIncreasing = Standard_True;
-              CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs, aPeriod, isIncreasing);
+              ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
+                                                      aPeriod, isIncreasing);
 
               //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
               //then after the next step (when U1 will be increased) U2 will be
@@ -2690,7 +2805,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
             }
           }
 
-          CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs, aV1[i], aV2[i]);
+          ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
+                                                              aV1[i], aV2[i]);
 
           if(isFirst)
           {
@@ -2771,11 +2887,9 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
             }
           }
 
-          AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
-                            theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
-                            anU1, aU2[i], aV1[i], aV1Prev[i],
-                            aV2[i], aV2Prev[i], i, isTheReverse,
-                            isForce, isFound1, isFound2);
+          aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i],
+                                      aV2[i], aV2Prev[i], i, isForce,
+                                      isFound1, isFound2);
 
           const Standard_Boolean isPrevVBound = !isVIntersect &&
                                                 ((Abs(aV1Prev[i] - aVSurf1f) <= theTol2D) ||
@@ -2783,7 +2897,6 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
                                                  (Abs(aV2Prev[i] - aVSurf2f) <= theTol2D) ||
                                                  (Abs(aV2Prev[i] - aVSurf2l) <= theTol2D));
 
-
           aV1Prev[i] = aV1[i];
           aV2Prev[i] = aV2[i];
 
@@ -2798,7 +2911,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
               aWLFindStatus[i] = WLFStatus_Exist;
             }
 
-            if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
+            if( (aWLFindStatus[i] != WLFStatus_Broken) ||
+                (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
             {
               if(aWLine[i]->NbPnts() > 0)
               {
@@ -2872,11 +2986,9 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
 
             Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
 
-            AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
-                              theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
-                              anU1, aU2[i], aV1[i], aV1Prev[i],
-                              aV2[i], aV2Prev[i], i, isTheReverse,
-                              Standard_False, isFound1, isFound2);
+            aBoundWork.AddBoundaryPoint(aWLine[i], anU1, aU2[i], aV1[i], aV1Prev[i],
+                                        aV2[i], aV2Prev[i], i,
+                                        Standard_False, isFound1, isFound2);
 
             if(isFound1 || isFound2)
             {
@@ -2952,7 +3064,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
                 continue;
               }
 
-              CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs, aU2[i], aV1[i], aV2[i]);
+              ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
+                                                                aU2[i], aV1[i], aV2[i]);
 
               AddPointIntoWL( theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
                               Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
@@ -2969,8 +3082,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
         //Step computing
 
         {
-          const Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
-                              aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
+          const Standard_Real aDeltaV1 = aRangeS1.Delta()/IntToReal(aNbPoints),
+                              aDeltaV2 = aRangeS2.Delta()/IntToReal(aNbPoints);
         
           math_Matrix aMatr(1, 3, 1, 5);
 
@@ -3078,8 +3191,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
             isAddedIntoWL[i] = Standard_True;
             SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), 
                                   anEquationCoeffs, i, aNbPoints, 1,
-                                  aWLine[i]->NbPnts(), aUSurf2f, aUSurf2l,
-                                  theTol2D, aPeriod, isTheReverse);
+                                  aWLine[i]->NbPnts(), theTol2D, aPeriod,
+                                  isTheReverse);
 
             aWLine[i]->ComputeVertexParameters(theTol3D);
             theSlin.Append(aWLine[i]);
@@ -3090,7 +3203,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
           isAddedIntoWL[i] = Standard_False;
         }
 
-#ifdef OCCT_DEBUG
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
         //aWLine[i]->Dump();
 #endif
       }
@@ -3104,7 +3217,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
   {
     for(Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
     {
-      Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::DownCast(theSlin.Value(aNbLin)));
+      Handle(IntPatch_WLine) aWLine1 (Handle(IntPatch_WLine)::
+                                            DownCast(theSlin.Value(aNbLin)));
 
       const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
       const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
@@ -3152,7 +3266,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
     for (Standard_Integer i = 0; i < aNbWLines; i++)
     {
       Standard_Real anU2t = 0.0;
-      if(!CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
+      if(!ComputationMethods::CylCylComputeParameters(anUl, i, anEquationCoeffs, anU2t))
         continue;
 
       const Standard_Real aDU2 = Abs(anU2t - aCurU2);
@@ -3168,7 +3282,7 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
     {
       Standard_Real anU2 = 0.0, anV1 = 0.0, anV2 = 0.0;
       Standard_Boolean isDone = 
-            CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
+            ComputationMethods::CylCylComputeParameters(anUf, anIndex, anEquationCoeffs,
                                     anU2, anV1, anV2);
       anUf += aDU;
 
@@ -3177,8 +3291,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
         continue;
       }
 
-      if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse, Standard_True,
-                        gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2),
+      if(AddPointIntoWL(theQuad1, theQuad2, anEquationCoeffs, isTheReverse,
+                        Standard_True, gp_Pnt2d(anUf, anV1), gp_Pnt2d(anU2, anV2),
                         aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
                         aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
                         aPeriod, aWLine->Curve(), anIndex, theTol3D,
@@ -3192,8 +3306,8 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
     {
       SeekAdditionalPoints( theQuad1, theQuad2, aWLine->Curve(), 
                             anEquationCoeffs, anIndex, aNbMinPoints,
-                            1, aWLine->NbPnts(), aUSurf2f, aUSurf2l,
-                            theTol2D, aPeriod, isTheReverse);
+                            1, aWLine->NbPnts(), theTol2D, aPeriod,
+                            isTheReverse);
 
       aWLine->ComputeVertexParameters(theTol3D);
       theSlin.Append(aWLine);
@@ -3614,7 +3728,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
            ptvalid = curvsol.Value(para);
            alig = new IntPatch_ALine(curvsol,Standard_False);
            kept = Standard_True;
-           //-- cout << "Transition indeterminee" << endl;
+           //-- std::cout << "Transition indeterminee" << std::endl;
          }
          if (kept) {
            Standard_Boolean Nfirstp = !firstp;
@@ -3696,46 +3810,3 @@ Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
   //
   return bFind;
 }
-//=======================================================================
-//function : IsToReverse
-//purpose  : 
-//=======================================================================
-Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
-                            const gp_Cylinder& Cy2,
-                            const Standard_Real Tol)
-{
-  Standard_Boolean bRet;
-  Standard_Real aR1,aR2, dR, aSc1, aSc2;
-  //
-  bRet=Standard_False;
-  //
-  aR1=Cy1.Radius();
-  aR2=Cy2.Radius();
-  dR=aR1-aR2;
-  if (dR<0.) {
-    dR=-dR;
-  }
-  if (dR>Tol) {
-    bRet=aR1>aR2;
-  }
-  else {
-    gp_Dir aDZ(0.,0.,1.);
-    //
-    const gp_Dir& aDir1=Cy1.Axis().Direction();
-    aSc1=aDir1*aDZ;
-    if (aSc1<0.) {
-      aSc1=-aSc1;
-    }
-    //
-    const gp_Dir& aDir2=Cy2.Axis().Direction();
-    aSc2=aDir2*aDZ;
-    if (aSc2<0.) {
-      aSc2=-aSc2;
-    }
-    //
-    if (aSc2<aSc1) {
-      bRet=!bRet;
-    }
-  }//if (dR<Tol) {
-  return bRet;
-}
index 75b6ba4..41d4dc0 100644 (file)
@@ -917,18 +917,9 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
     ListOfPnts.Clear();
     if(isGeomInt)
     {
-      if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite())
-      {
-        GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, 
-                        TolTang, ListOfPnts, RestrictLine,
-                        typs1, typs2, theIsReqToKeepRLine);
-      }
-      else
-      {
-        GeomGeomPerfomTrimSurf( theS1, theD1, theS2, theD2,
-                                TolArc, TolTang, ListOfPnts, RestrictLine,
-                                typs1, typs2, theIsReqToKeepRLine);
-      }
+      GeomGeomPerfom( theS1, theD1, theS2, theD2, TolArc, 
+                      TolTang, ListOfPnts, RestrictLine,
+                      typs1, typs2, theIsReqToKeepRLine);
     }
     else
     {
@@ -1182,16 +1173,8 @@ void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  theS1,
   }
   else if(ts1 == 1)
   {
-    if(theD1->DomainIsInfinite() || theD2->DomainIsInfinite())
-    {
-      GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, 
-                      TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine);
-    }
-    else
-    {
-      GeomGeomPerfomTrimSurf(theS1, theD1, theS2, theD2,
-              TolArc, TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine);
-    }
+    GeomGeomPerfom(theS1, theD1, theS2, theD2, TolArc, 
+                    TolTang, ListOfPnts, RestrictLine, typs1, typs2, theIsReqToKeepRLine);
   }
 
   if(!theIsReqToPostWLProc)
@@ -1438,7 +1421,36 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the
           slin.Append(wlin);
         }
         else
+        {
+          if(line->ArcType() == IntPatch_Walking)
+          {
+            Handle(IntPatch_WLine)::DownCast(line)->EnablePurging(Standard_False);
+          }
+
           slin.Append(line);
+        }
+      }
+
+      for (Standard_Integer i = 1; i <= interii.NbPnts(); i++)
+      {
+        spnt.Append(interii.Point(i));
+      }
+
+      if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder))
+      {
+        IntPatch_WLineTool::JoinWLines( slin, spnt, TolTang,
+                                        theS1->IsUPeriodic()? theS1->UPeriod() : 0.0,
+                                        theS2->IsUPeriodic()? theS2->UPeriod() : 0.0,
+                                        theS1->IsVPeriodic()? theS1->VPeriod() : 0.0,
+                                        theS2->IsVPeriodic()? theS2->VPeriod() : 0.0,
+                                        theS1->FirstUParameter(),
+                                        theS1->LastUParameter(),
+                                        theS1->FirstVParameter(),
+                                        theS1->LastVParameter(),
+                                        theS2->FirstUParameter(),
+                                        theS2->LastUParameter(),
+                                        theS2->FirstVParameter(),
+                                        theS2->LastVParameter());
       }
 
       if(isQuadSet)
@@ -1472,11 +1484,6 @@ void IntPatch_Intersection::GeomGeomPerfom(const Handle(Adaptor3d_HSurface)& the
                                      theS2->IsVPeriodic()? theS2->VPeriod() : 0.0,
                                      aBx1, aBx2);
       }
-
-      for (Standard_Integer i = 1; i <= interii.NbPnts(); i++)
-      {
-        spnt.Append(interii.Point(i));
-      }
     }
   }
   else
@@ -1568,84 +1575,6 @@ void IntPatch_Intersection::
   }
 }
 
-//=======================================================================
-//function : GeomGeomPerfomTrimSurf
-//purpose  : This function returns ready walking-line (which is not need
-//            in convertation) as an intersection line between two
-//            trimmed surfaces.
-//=======================================================================
-void IntPatch_Intersection::
-  GeomGeomPerfomTrimSurf( const Handle(Adaptor3d_HSurface)& theS1,
-                          const Handle(Adaptor3d_TopolTool)& theD1,
-                          const Handle(Adaptor3d_HSurface)& theS2,
-                          const Handle(Adaptor3d_TopolTool)& theD2,
-                          const Standard_Real theTolArc,
-                          const Standard_Real theTolTang,
-                          IntSurf_ListOfPntOn2S& theListOfPnts,
-                          const Standard_Boolean RestrictLine,
-                          const GeomAbs_SurfaceType theTyps1,
-                          const GeomAbs_SurfaceType theTyps2,
-                          const Standard_Boolean theIsReqToKeepRLine)
-{
-  IntSurf_Quadric Quad1,Quad2;
-
-  if((theTyps1 == GeomAbs_Cylinder) && (theTyps2 == GeomAbs_Cylinder))
-  {
-    IntPatch_ImpImpIntersection anInt;
-    anInt.Perform(theS1, theD1, theS2, theD2, myTolArc,
-                  myTolTang, Standard_True, theIsReqToKeepRLine);
-
-    done = anInt.IsDone();
-
-    if(done)
-    {
-      empt = anInt.IsEmpty();
-      if (!empt)
-      {
-        tgte = anInt.TangentFaces();
-        if (tgte)
-          oppo = anInt.OppositeFaces();
-
-        const Standard_Integer aNbLin = anInt.NbLines();
-        const Standard_Integer aNbPts = anInt.NbPnts();
-
-        for(Standard_Integer aLID = 1; aLID <= aNbLin; aLID++)
-        {
-          const Handle(IntPatch_Line)& aLine = anInt.Line(aLID);
-          slin.Append(aLine);
-        }
-
-        for(Standard_Integer aPID = 1; aPID <= aNbPts; aPID++)
-        {
-          const IntPatch_Point& aPoint = anInt.Point(aPID);
-          spnt.Append(aPoint);
-        }
-
-        IntPatch_WLineTool::JoinWLines( slin, spnt, theTolTang,
-                                        theS1->IsUPeriodic()? theS1->UPeriod() : 0.0,
-                                        theS2->IsUPeriodic()? theS2->UPeriod() : 0.0,
-                                        theS1->IsVPeriodic()? theS1->VPeriod() : 0.0,
-                                        theS2->IsVPeriodic()? theS2->VPeriod() : 0.0,
-                                        theS1->FirstUParameter(),
-                                        theS1->LastUParameter(),
-                                        theS1->FirstVParameter(),
-                                        theS1->LastVParameter(),
-                                        theS2->FirstUParameter(),
-                                        theS2->LastUParameter(),
-                                        theS2->FirstVParameter(),
-                                        theS2->LastVParameter());
-      }
-    }
-  }
-  else
-  {
-    GeomGeomPerfom(theS1, theD1, theS2, theD2,
-            theTolArc, theTolTang, theListOfPnts,
-            RestrictLine, theTyps1, theTyps2, theIsReqToKeepRLine);
-  }
-}
-
-
 void IntPatch_Intersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
                                     const Handle(Adaptor3d_TopolTool)& D1,
                                     const Handle(Adaptor3d_HSurface)&  S2,
index cfe6b03..96a4dac 100644 (file)
@@ -167,14 +167,6 @@ private:
   //! will always keep both lines even if they are coincided.
   Standard_EXPORT void GeomGeomPerfom (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, IntSurf_ListOfPntOn2S& LOfPnts, const Standard_Boolean RestrictLine, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2, const Standard_Boolean theIsReqToKeepRLine = Standard_False);
   
-  //! Flag theIsReqToKeepRLine has been enterred only for
-  //! compatibility with TopOpeBRep package. It shall be deleted
-  //! after deleting TopOpeBRep.
-  //! When intersection result returns IntPatch_RLine and another
-  //! IntPatch_Line (not restriction) we (in case of theIsReqToKeepRLine==TRUE)
-  //! will always keep both lines even if they are coincided.
-  Standard_EXPORT void GeomGeomPerfomTrimSurf (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Real TolArc, const Standard_Real TolTang, IntSurf_ListOfPntOn2S& LOfPnts, const Standard_Boolean RestrictLine, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2, const Standard_Boolean theIsReqToKeepRLine = Standard_False);
-  
   Standard_EXPORT void GeomParamPerfom (const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_TopolTool)& D1, const Handle(Adaptor3d_HSurface)& S2, const Handle(Adaptor3d_TopolTool)& D2, const Standard_Boolean isNotAnalitical, const GeomAbs_SurfaceType typs1, const GeomAbs_SurfaceType typs2);
 
 
index ca38c1c..bb09d22 100644 (file)
@@ -17,6 +17,9 @@
 #include <Adaptor3d_TopolTool.hxx>
 #include <ElCLib.hxx>
 
+// It is pure empirical value.
+const Standard_Real IntPatch_WLineTool::myMaxConcatAngle = M_PI/6;
+
 //Bit-mask is used for information about 
 //the operation made in
 //IntPatch_WLineTool::ExtendTwoWlinesToEachOther() method.
@@ -830,27 +833,26 @@ static Standard_Boolean IsOutOfDomain(const Bnd_Box2d& theBoxS1,
 }
 
 //=======================================================================
-//function : CheckArguments
+//function : CheckArgumentsToExtend
 //purpose  : Check if extending is possible
 //            (see IntPatch_WLineTool::ExtendTwoWlinesToEachOther)
 //=======================================================================
-static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1,
-                                       const IntSurf_Quadric& theS2,
-                                       const IntSurf_PntOn2S& thePtWL1,
-                                       const IntSurf_PntOn2S& thePtWL2,
-                                       IntSurf_PntOn2S& theNewPoint,
-                                       const gp_Vec& theVec1,
-                                       const gp_Vec& theVec2,
-                                       const gp_Vec& theVec3,
-                                       const Bnd_Box2d& theBoxS1,
-                                       const Bnd_Box2d& theBoxS2,
-                                       const Standard_Real theToler3D,
-                                       const Standard_Real theU1Period,
-                                       const Standard_Real theU2Period,
-                                       const Standard_Real theV1Period,
-                                       const Standard_Real theV2Period)
+Standard_Boolean CheckArgumentsToExtend(const IntSurf_Quadric& theS1,
+                                        const IntSurf_Quadric& theS2,
+                                        const IntSurf_PntOn2S& thePtWL1,
+                                        const IntSurf_PntOn2S& thePtWL2,
+                                        IntSurf_PntOn2S& theNewPoint,
+                                        const gp_Vec& theVec1,
+                                        const gp_Vec& theVec2,
+                                        const gp_Vec& theVec3,
+                                        const Bnd_Box2d& theBoxS1,
+                                        const Bnd_Box2d& theBoxS2,
+                                        const Standard_Real theToler3D,
+                                        const Standard_Real theU1Period,
+                                        const Standard_Real theU2Period,
+                                        const Standard_Real theV1Period,
+                                        const Standard_Real theV2Period)
 {
-  const Standard_Real aMaxAngle = M_PI/6; //30 degree
   const Standard_Real aSqToler = theToler3D*theToler3D;
 
   if(theVec3.SquareMagnitude() <= aSqToler)
@@ -858,9 +860,9 @@ static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1,
     return Standard_False;
   }
 
-  if((theVec1.Angle(theVec2) > aMaxAngle) ||
-     (theVec1.Angle(theVec3) > aMaxAngle) ||
-     (theVec2.Angle(theVec3) > aMaxAngle))
+  if((theVec1.Angle(theVec2) > IntPatch_WLineTool::myMaxConcatAngle) ||
+     (theVec1.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle) ||
+     (theVec2.Angle(theVec3) > IntPatch_WLineTool::myMaxConcatAngle))
   {
     return Standard_False;
   }
@@ -908,6 +910,20 @@ static Standard_Boolean CheckArguments(const IntSurf_Quadric& theS1,
 }
 
 //=======================================================================
+//function : CheckArgumentsToJoin
+//purpose  : Check if joining is possible
+//            (see IntPatch_WLineTool::JoinWLines)
+//=======================================================================
+Standard_Boolean CheckArgumentsToJoin(const gp_Vec& theVec1,
+                                      const gp_Vec& theVec2)
+{
+  // [0, PI] - range
+  const Standard_Real anAngle = theVec1.Angle(theVec2);
+
+  return (anAngle < IntPatch_WLineTool::myMaxConcatAngle);
+}
+
+//=======================================================================
 //function : ExtendTwoWLFirstFirst
 //purpose  : Performs extending theWLine1 and theWLine2 through their
 //            respecting start point.
@@ -932,10 +948,10 @@ static void ExtendTwoWLFirstFirst(const IntSurf_Quadric& theS1,
                                    Standard_Boolean &theHasBeenJoined)
 {
   IntSurf_PntOn2S aPOn2S;
-  if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
-                     theVec1, theVec2, theVec3, 
-                     theBoxS1, theBoxS2, theToler3D,
-                     theU1Period, theU2Period, theV1Period, theV2Period))
+  if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
+                             theVec1, theVec2, theVec3, 
+                             theBoxS1, theBoxS2, theToler3D,
+                             theU1Period, theU2Period, theV1Period, theV2Period))
   {
     return;
   }
@@ -1006,10 +1022,10 @@ static void ExtendTwoWLFirstLast(const IntSurf_Quadric& theS1,
                                  Standard_Boolean &theHasBeenJoined)
 {
   IntSurf_PntOn2S aPOn2S;
-  if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
-                     theVec1, theVec2, theVec3, 
-                     theBoxS1, theBoxS2, theToler3D,
-                     theU1Period, theU2Period, theV1Period, theV2Period))
+  if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
+                             theVec1, theVec2, theVec3, 
+                             theBoxS1, theBoxS2, theToler3D,
+                             theU1Period, theU2Period, theV1Period, theV2Period))
   {
     return;
   }
@@ -1078,10 +1094,10 @@ static void ExtendTwoWLLastFirst(const IntSurf_Quadric& theS1,
                                  Standard_Boolean &theHasBeenJoined)
 {
   IntSurf_PntOn2S aPOn2S;
-  if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
-                     theVec1, theVec2, theVec3, 
-                     theBoxS1, theBoxS2, theToler3D,
-                     theU1Period, theU2Period, theV1Period, theV2Period))
+  if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
+                             theVec1, theVec2, theVec3, 
+                             theBoxS1, theBoxS2, theToler3D,
+                             theU1Period, theU2Period, theV1Period, theV2Period))
   {
     return;
   }
@@ -1146,10 +1162,10 @@ static void ExtendTwoWLLastLast(const IntSurf_Quadric& theS1,
                                 Standard_Boolean &theHasBeenJoined)
 {
   IntSurf_PntOn2S aPOn2S;
-  if(!CheckArguments(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
-                     theVec1, theVec2, theVec3, 
-                     theBoxS1, theBoxS2, theToler3D,
-                     theU1Period, theU2Period, theV1Period, theV2Period))
+  if(!CheckArgumentsToExtend(theS1, theS2, thePtWL1, thePtWL2, aPOn2S,
+                             theVec1, theVec2, theVec3, 
+                             theBoxS1, theBoxS2, theToler3D,
+                             theU1Period, theU2Period, theV1Period, theV2Period))
   {
     return;
   }
@@ -1387,10 +1403,19 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       {
         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2);
-        if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period,
-                          theV1Period, theV2Period, theUfSurf1, theUlSurf1,
-                          theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
-                          theVfSurf2, theVlSurf2))
+        Standard_Boolean aCond = 
+              CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()),
+                                   gp_Vec(aPt2.Value(), aPntFWL2.Value()));
+
+        aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1,
+                                        theU1Period, theU2Period,
+                                        theV1Period, theV2Period,
+                                        theUfSurf1, theUlSurf1,
+                                        theVfSurf1, theVlSurf1,
+                                        theUfSurf2, theUlSurf2,
+                                        theVfSurf2, theVlSurf2);
+          
+        if(aCond)
         {
           aWLine1->ClearVertexes();
           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
@@ -1413,10 +1438,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       {
         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(2);
         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1);
-        if(!IsSeamOrBound(aPt1, aPt2, aPntFWL1, theU1Period, theU2Period,
-                          theV1Period, theV2Period, theUfSurf1, theUlSurf1,
-                          theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
-                          theVfSurf2, theVlSurf2))
+
+        Standard_Boolean aCond = 
+              CheckArgumentsToJoin(gp_Vec(aPntFWL1.Value(), aPt1.Value()),
+                                   gp_Vec(aPt2.Value(), aPntLWL2.Value()));
+
+        aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntFWL1,
+                                        theU1Period, theU2Period,
+                                        theV1Period, theV2Period,
+                                        theUfSurf1, theUlSurf1,
+                                        theVfSurf1, theVlSurf1,
+                                        theUfSurf2, theUlSurf2,
+                                        theVfSurf2, theVlSurf2);
+
+        if(aCond)
         {
           aWLine1->ClearVertexes();
           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
@@ -1439,10 +1474,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       {
         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1);
         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(2);
-        if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period,
-                          theV1Period, theV2Period, theUfSurf1, theUlSurf1,
-                          theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
-                          theVfSurf2, theVlSurf2))
+
+        Standard_Boolean aCond = 
+              CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()),
+                                   gp_Vec(aPntFWL2.Value(), aPt2.Value()));
+
+        aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1,
+                                        theU1Period, theU2Period,
+                                        theV1Period, theV2Period,
+                                        theUfSurf1, theUlSurf1,
+                                        theVfSurf1, theVlSurf1,
+                                        theUfSurf2, theUlSurf2,
+                                        theVfSurf2, theVlSurf2);
+
+        if(aCond)
         {
           aWLine1->ClearVertexes();
           for(Standard_Integer aNPt = 1; aNPt <= aNbPntsWL2; aNPt++)
@@ -1465,10 +1510,20 @@ void IntPatch_WLineTool::JoinWLines(IntPatch_SequenceOfLine& theSlin,
       {
         const IntSurf_PntOn2S& aPt1 = aWLine1->Point(aNbPntsWL1-1);
         const IntSurf_PntOn2S& aPt2 = aWLine2->Point(aNbPntsWL2-1);
-        if(!IsSeamOrBound(aPt1, aPt2, aPntLWL1, theU1Period, theU2Period,
-                          theV1Period, theV2Period, theUfSurf1, theUlSurf1,
-                          theVfSurf1, theVlSurf1, theUfSurf2, theUlSurf2,
-                          theVfSurf2, theVlSurf2))
+
+        Standard_Boolean aCond = 
+              CheckArgumentsToJoin(gp_Vec(aPt1.Value(), aPntLWL1.Value()),
+                                   gp_Vec(aPntLWL2.Value(), aPt2.Value()));
+
+        aCond = aCond && !IsSeamOrBound(aPt1, aPt2, aPntLWL1, 
+                                        theU1Period, theU2Period,
+                                        theV1Period, theV2Period,
+                                        theUfSurf1, theUlSurf1,
+                                        theVfSurf1, theVlSurf1,
+                                        theUfSurf2, theUlSurf2,
+                                        theVfSurf2, theVlSurf2);
+
+        if(aCond)
         {
           aWLine1->ClearVertexes();
           for(Standard_Integer aNPt = aNbPntsWL2; aNPt >= 1; aNPt--)
index 90b0d5d..02e4584 100644 (file)
@@ -96,6 +96,9 @@ public:
                                                          const Standard_Real theV2Period,
                                                          const Bnd_Box2d& theBoxS1,
                                                          const Bnd_Box2d& theBoxS2);
+
+//! Max angle to concatenate two WLines to avoid result with C0-continuity
+  static const Standard_Real myMaxConcatAngle;
 };
 
 #endif
\ No newline at end of file
index 1f3e4ad..1c2a537 100644 (file)
@@ -3,4 +3,4 @@ tscale s 0 0 0 SCALE
 explode s E
 blend result s SCALE*2 s_5
  
-checkprops result -s 1.65391e+08 -eps 0.1
+checkprops result -s 1.65391e+008
index 7e4b81e..6dad90a 100755 (executable)
@@ -31,16 +31,19 @@ explode b1 f
 explode b2 f
 
 smallview
-donly b1_4 b2_1
+donly b1_4
 fit
+display b2_1
 
 dchrono h reset
 dchrono h start
 
 bopcurves b1_4 b2_1 -2d
 
 dchrono h stop
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_1.png
+
 set q [dchrono h show]
 
 regexp {CPU user time: ([-0-9.+eE]+) seconds} $q full z
@@ -56,10 +59,10 @@ if { $z > ${max_time} } {
 mksurface s1 b1_4
 mksurface s2 b2_1
 
-dchrono h2 stop
-set q2 [dchrono h2 show]
+dchrono h2 reset
+dchrono h2 start
 
-set CurveNumb [intersect i s1 s2]
+set CurveNumb [intersect resi s1 s2]
 
 dchrono h2 stop
 set q2 [dchrono h2 show]
@@ -79,4 +82,8 @@ if { [llength ${CurveNumb}] < 1 } {
     puts "OK : Good intersection"
 }
 
-checkview -screenshot -2d -path ${imagedir}/${test_image}.png
+don resi*
+fit
+display s1 s2
+
+checkview -screenshot -2d -path ${imagedir}/${test_image}_2.png
diff --git a/tests/bugs/modalg_6/bug23178 b/tests/bugs/modalg_6/bug23178
new file mode 100644 (file)
index 0000000..a16ee98
--- /dev/null
@@ -0,0 +1,43 @@
+puts "================"
+puts "OCC23178"
+puts "================"
+puts ""
+#######################################################################
+#  Intersection of cylinders fails to produce results
+#######################################################################
+
+set GoodNbCurv 6
+
+foreach c [directory result*] {
+  unset $c
+}
+
+cylinder s1 4.999999999813, 1.35836143386791, 0.20975890778078 -1, 0, 0 0, 0, 1 10
+cylinder s2 9.999999999813, 9.1401960568287, 0.11837300583107 0, -0.75870713028692, 0.651431877061437 0, -0.651431877061437, -0.75870713028692 5 
+
+intersect result s1 s2
+
+foreach c [directory result*] {
+  bounds $c U1 U2
+  
+  if {[dval U2-U1] < 1.0e-9} {
+    puts "Error: Wrong curve's range!"
+  }
+  
+  xdistcs $c s1 U1 U2 10 2.0e-7
+  xdistcs $c s2 U1 U2 10 2.0e-7    
+}
+
+set NbCurv [llength [directory result*]]
+
+if { $NbCurv == $GoodNbCurv } {
+  puts "OK: Number of curves is good!"
+} else {
+  puts "Error: $GoodNbCurv is expected but $NbCurv is found!"
+}
+
+smallview
+don result*
+fit
+disp s1 s2
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png
index fe5faee..8c57872 100644 (file)
@@ -6,7 +6,7 @@ puts ""
 # Huge tolerance obtained in the result of intersection of two cylindrical faces
 #################################################
 
-set ExpTol 0.00027580830315682321
+set ExpTol 6.9230965382090094e-006
 set GoodNbCurv 2
 
 restore [locate_data_file OCC496a.brep] a 
index d992b88..c039cd8 100644 (file)
@@ -7,7 +7,7 @@ puts ""
 #################################################
 
 set ExpTol 1.0e-7
-set GoodNbCurv 5
+set GoodNbCurv 3
 
 restore [locate_data_file bug27761_c1.brep] c1
 restore [locate_data_file bug27761_c2.brep] c2