0028599: Replacement of old Boolean operations with new ones in BRepProj_Projection...
[occt.git] / src / IntPolyh / IntPolyh_MaillageAffinage.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 4affa9d..a0f7958
@@ -1,52 +1,58 @@
 // Created on: 1999-03-05
 // Created by: Fabrice SERVANT
 // Copyright (c) 1999-1999 Matra Datavision
-// Copyright (c) 1999-2012 OPEN CASCADE SAS
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
 //
-// The content of this file is subject to the Open CASCADE Technology Public
-// License Version 6.5 (the "License"). You may not use the content of this file
-// except in compliance with the License. Please obtain a copy of the License
-// at http://www.opencascade.org and read it completely before using this file.
+// This file is part of Open CASCADE Technology software library.
 //
-// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
-// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
+// 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.
 //
-// The Original Code and all software distributed under the License is
-// distributed on an "AS IS" basis, without warranty of any kind, and the
-// Initial Developer hereby disclaims all such warranties, including without
-// limitation, any warranties of merchantability, fitness for a particular
-// purpose or non-infringement. Please see the License for the specific terms
-// and conditions governing the rights and limitations under the License.
-
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
 
 //  modified by Edward AGAPOV (eap) Tue Jan 22 2002 (bug occ53)
 //  - improve SectionLine table management (avoid memory reallocation)
 //  - some protection against arrays overflow
-
 //  modified by Edward AGAPOV (eap) Thu Feb 14 2002 (occ139)
 //  - make Section Line parts rightly connected (prepend 2nd part to the 1st)
 //  - TriangleShape() for debugging purpose
-
 //  Modified by skv - Thu Sep 25 17:42:42 2003 OCC567
 //  modified by ofv Thu Apr  8 14:58:13 2004 fip
 
-
-#include <IntPolyh_MaillageAffinage.ixx>
-
-#include <Precision.hxx>
-#include <gp_Pnt.hxx>
+#include <Adaptor3d_HSurface.hxx>
+#include <Bnd_BoundSortBox.hxx>
+#include <Bnd_Box.hxx>
+#include <Bnd_HArray1OfBox.hxx>
 #include <gp.hxx>
-
+#include <gp_Pnt.hxx>
+#include <IntPolyh_ListOfCouples.hxx>
+#include <IntPolyh_Couple.hxx>
+#include <IntPolyh_Edge.hxx>
+#include <IntPolyh_MaillageAffinage.hxx>
+#include <IntPolyh_Point.hxx>
+#include <IntPolyh_SectionLine.hxx>
+#include <IntPolyh_StartPoint.hxx>
+#include <IntPolyh_Tools.hxx>
+#include <IntPolyh_Triangle.hxx>
+#include <Precision.hxx>
+#include <TColStd_Array1OfInteger.hxx>
+#include <TColStd_MapOfInteger.hxx>
 #include <TColStd_ListIteratorOfListOfInteger.hxx>
+#include <NCollection_UBTree.hxx>
+#include <NCollection_UBTreeFiller.hxx>
+#include <algorithm>
+#include <NCollection_IndexedDataMap.hxx>
 
-#include <Bnd_BoundSortBox.hxx>
-#include <Bnd_HArray1OfBox.hxx> 
-
-#include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
+typedef NCollection_Array1<Standard_Integer> IntPolyh_ArrayOfInteger;
+typedef NCollection_IndexedDataMap
+  <Standard_Integer,
+   IntPolyh_ArrayOfInteger,
+   TColStd_MapIntegerHasher> IntPolyh_IndexedDataMapOfIntegerArrayOfInteger;
 
-#include <IntPolyh_ArrayOfCouples.hxx>
-#include <IntPolyh_Edge.hxx>
-#include <IntPolyh_Couple.hxx>
 
 static Standard_Real MyTolerance=10.0e-7;
 static Standard_Real MyConfusionPrecision=10.0e-12;
@@ -54,92 +60,171 @@ static Standard_Real SquareMyConfusionPrecision=10.0e-24;
 //
 static
   inline Standard_Real maxSR(const Standard_Real a,
-                            const Standard_Real b,
-                            const Standard_Real c);
+                             const Standard_Real b,
+                             const Standard_Real c);
 
 static
   inline Standard_Real minSR(const Standard_Real a,
-                            const Standard_Real b,
-                            const Standard_Real c);
+                             const Standard_Real b,
+                             const Standard_Real c);
 static
   Standard_Integer project6(const IntPolyh_Point &ax, 
-                           const IntPolyh_Point &p1,
-                           const IntPolyh_Point &p2, 
-                           const IntPolyh_Point &p3, 
-                           const IntPolyh_Point &q1,
-                           const IntPolyh_Point &q2, 
-                           const IntPolyh_Point &q3);
+                            const IntPolyh_Point &p1,
+                            const IntPolyh_Point &p2, 
+                            const IntPolyh_Point &p3, 
+                            const IntPolyh_Point &q1,
+                            const IntPolyh_Point &q2, 
+                            const IntPolyh_Point &q3);
 static
   void TestNbPoints(const Standard_Integer ,
-                   Standard_Integer &NbPoints,
-                   Standard_Integer &NbPointsTotal,
-                   const IntPolyh_StartPoint &Pt1,
-                   const IntPolyh_StartPoint &Pt2,
-                   IntPolyh_StartPoint &SP1,
-                   IntPolyh_StartPoint &SP2);
+                    Standard_Integer &NbPoints,
+                    Standard_Integer &NbPointsTotal,
+                    const IntPolyh_StartPoint &Pt1,
+                    const IntPolyh_StartPoint &Pt2,
+                    IntPolyh_StartPoint &SP1,
+                    IntPolyh_StartPoint &SP2);
 static
   void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
-                                       const IntPolyh_Point &NormaleTri,
-                                       const IntPolyh_Point &PE1,
-                                       const IntPolyh_Point &PE2,
-                                       const IntPolyh_Point &Edge,
-                                       const IntPolyh_Point &PT1,
-                                       const IntPolyh_Point &PT2,
-                                       const IntPolyh_Point &Cote,
-                                       const Standard_Integer CoteIndex,
-                                       IntPolyh_StartPoint &SP1,
-                                       IntPolyh_StartPoint &SP2,
-                                       Standard_Integer &NbPoints);
-static
-  void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
-                                        const IntPolyh_Point &NormaleTri,
-                                        const IntPolyh_Triangle &Tri1,
-                                        const IntPolyh_Triangle &Tri2,
-                                        const IntPolyh_Point &PE1,
-                                        const IntPolyh_Point &PE2,
-                                        const IntPolyh_Point &Edge,
-                                        const Standard_Integer EdgeIndex,
-                                        const IntPolyh_Point &PT1,
-                                        const IntPolyh_Point &PT2,
-                                        const IntPolyh_Point &Cote,
-                                      const Standard_Integer CoteIndex,
-                                        IntPolyh_StartPoint &SP1,
-                                        IntPolyh_StartPoint &SP2,
-                                        Standard_Integer &NbPoints); 
+                                        const IntPolyh_Point &NormaleTri,
+                                        const IntPolyh_Triangle &Tri1,
+                                        const IntPolyh_Triangle &Tri2,
+                                        const IntPolyh_Point &PE1,
+                                        const IntPolyh_Point &PE2,
+                                        const IntPolyh_Point &Edge,
+                                        const Standard_Integer EdgeIndex,
+                                        const IntPolyh_Point &PT1,
+                                        const IntPolyh_Point &PT2,
+                                        const IntPolyh_Point &Cote,
+                                        const Standard_Integer CoteIndex,
+                                        IntPolyh_StartPoint &SP1,
+                                        IntPolyh_StartPoint &SP2,
+                                        Standard_Integer &NbPoints); 
 static
-  Standard_Integer CheckCoupleAndGetAngle(const Standard_Integer T1, 
-                                         const Standard_Integer T2,
-                                         Standard_Real& Angle, 
-                                         IntPolyh_ArrayOfCouples &TTrianglesContacts);
+  Standard_Boolean CheckCoupleAndGetAngle(const Standard_Integer T1, 
+                                          const Standard_Integer T2,
+                                          Standard_Real& Angle, 
+                                          IntPolyh_ListOfCouples &TTrianglesContacts);
 static
-  Standard_Integer CheckCoupleAndGetAngle2(const Standard_Integer T1,
-                                          const Standard_Integer T2,
-                                          const Standard_Integer T11, 
-                                          const Standard_Integer T22,
-                                          Standard_Integer &CT11,
-                                          Standard_Integer &CT22, 
-                                          Standard_Real & Angle,
-                                          IntPolyh_ArrayOfCouples &TTrianglesContacts);
+  Standard_Boolean CheckCoupleAndGetAngle2(const Standard_Integer T1,
+                                           const Standard_Integer T2,
+                                           const Standard_Integer T11, 
+                                           const Standard_Integer T22,
+                                           IntPolyh_ListIteratorOfListOfCouples& theItCT11,
+                                           IntPolyh_ListIteratorOfListOfCouples& theItCT22,
+                                           Standard_Real & Angle,
+                                           IntPolyh_ListOfCouples &TTrianglesContacts);
 static
   Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
-                                      IntPolyh_ArrayOfTangentZones & TTangentZones,
-                                      IntPolyh_StartPoint & SP,
-                                      const Standard_Boolean Prepend=Standard_False); 
+                                       IntPolyh_ArrayOfTangentZones & TTangentZones,
+                                       IntPolyh_StartPoint & SP,
+                                       const Standard_Boolean Prepend=Standard_False); 
 
-//modified by NIZNHY-PKV Fri Jan 20 11:01:30 2012f
 static
   Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
-                                const Standard_Integer aIndex,
-                                const Standard_Real aTol2,
-                                Standard_Real& aDegX);
+                                 const Standard_Integer aIndex,
+                                 const Standard_Real aTol2,
+                                 Standard_Real& aDegX);
 static
   void DegeneratedIndex(const TColStd_Array1OfReal& Xpars,
-                       const Standard_Integer aNbX,
-                       const Handle(Adaptor3d_HSurface)& aS,
-                       const Standard_Integer aIsoDirection,
-                       Standard_Integer& aI1,
-                       Standard_Integer& aI2);
-//modified by NIZNHY-PKV Fri Jan 20 11:01:32 2012t
+                        const Standard_Integer aNbX,
+                        const Handle(Adaptor3d_HSurface)& aS,
+                        const Standard_Integer aIsoDirection,
+                        Standard_Integer& aI1,
+                        Standard_Integer& aI2);
+
+//=======================================================================
+//class : IntPolyh_BoxBndTreeSelector
+//purpose  : Select interfering bounding boxes
+//=======================================================================
+typedef NCollection_UBTree<Standard_Integer, Bnd_Box> IntPolyh_BoxBndTree;
+class IntPolyh_BoxBndTreeSelector : public IntPolyh_BoxBndTree::Selector {
+ public:
+  IntPolyh_BoxBndTreeSelector(const Bnd_Box& theBox) : myBox(theBox) {}
+  //
+  virtual Standard_Boolean Reject(const Bnd_Box& theOther) const
+  {
+    return myBox.IsOut(theOther);
+  }
+  //
+  virtual Standard_Boolean Accept(const Standard_Integer &theInd)
+  {
+    myIndices.Append(theInd);
+    return Standard_True;
+  }
+  //
+  const TColStd_ListOfInteger& Indices() const
+  {
+    return myIndices;
+  }
+ private:
+  Bnd_Box myBox;
+  TColStd_ListOfInteger myIndices;
+};
+
+//=======================================================================
+//function : GetInterferingTriangles
+//purpose  : Returns indices of the triangles with interfering bounding boxes
+//=======================================================================
+static
+  void GetInterferingTriangles(IntPolyh_ArrayOfTriangles& theTriangles1,
+                               const IntPolyh_ArrayOfPoints& thePoints1,
+                               IntPolyh_ArrayOfTriangles& theTriangles2,
+                               const IntPolyh_ArrayOfPoints& thePoints2,
+                               IntPolyh_IndexedDataMapOfIntegerArrayOfInteger& theCouples)
+{
+  // To find the triangles with interfering bounding boxes
+  // use the algorithm of unbalanced binary tree of overlapping bounding boxes
+  IntPolyh_BoxBndTree aBBTree;
+  NCollection_UBTreeFiller <Standard_Integer, Bnd_Box> aTreeFiller(aBBTree);
+  // 1. Fill the tree with the boxes of the triangles from second surface
+  Standard_Integer i, aNbT2 = theTriangles2.NbItems();
+  Standard_Boolean bAdded = Standard_False;
+  for (i = 0; i < aNbT2; ++i) {
+    IntPolyh_Triangle& aT = theTriangles2[i];
+    if (!aT.IsIntersectionPossible() || aT.IsDegenerated()) {
+      continue;
+    }
+    //
+    const Bnd_Box& aBox = aT.BoundingBox(thePoints2);
+    aTreeFiller.Add(i, aBox);
+    bAdded = Standard_True;
+  }
+  //
+  if (!bAdded)
+    // Intersection is not possible for all triangles in theTriangles2
+    return;
+
+  // 2. Shake the tree filler
+  aTreeFiller.Fill();
+  //
+  // 3. Find boxes interfering with the first triangles
+  Standard_Integer aNbT1 = theTriangles1.NbItems();
+  for (i = 0; i < aNbT1; ++i) {
+    IntPolyh_Triangle& aT = theTriangles1[i];
+    if (!aT.IsIntersectionPossible() || aT.IsDegenerated()) {
+      continue;
+    }
+    //
+    const Bnd_Box& aBox = aT.BoundingBox(thePoints1);
+    //
+    IntPolyh_BoxBndTreeSelector aSelector(aBox);
+    if (!aBBTree.Select(aSelector)) {
+      continue;
+    }
+    //
+    const TColStd_ListOfInteger& aLI = aSelector.Indices();
+    // Sort the indices
+    IntPolyh_ArrayOfInteger anArr(1, aLI.Extent());
+    TColStd_ListIteratorOfListOfInteger aItLI(aLI);
+    for (Standard_Integer j = 1; aItLI.More(); aItLI.Next(), ++j) {
+      anArr(j) = aItLI.Value();
+    }
+    //
+    std::sort(anArr.begin(), anArr.end());
+    //
+    theCouples.Add(i, anArr);
+  }
+}
 
 //=======================================================================
 //function : IntPolyh_MaillageAffinage
@@ -160,19 +245,8 @@ IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
   FlecheMax2(0.0), 
   FlecheMin1(0.0), 
   FlecheMin2(0.0),
-  FlecheMoy1(0.0), 
-  FlecheMoy2(0.0), 
   myEnlargeZone(Standard_False) 
 { 
-   TPoints1.Init(10000);
-   TEdges1.Init(30000);
-   TTriangles1.Init(20000);
-   
-   TPoints2.Init(10000);
-   TEdges2.Init(30000);
-   TTriangles2.Init(20000);
-  
-   TStartPoints.Init(10000);
 }
 //=======================================================================
 //function : IntPolyh_MaillageAffinage
@@ -197,20 +271,23 @@ IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
   FlecheMax2(0.0), 
   FlecheMin1(0.0), 
   FlecheMin2(0.0),
-  FlecheMoy1(0.0), 
-  FlecheMoy2(0.0), 
   myEnlargeZone(Standard_False)
 { 
-   TPoints1.Init(10000);
-   TEdges1.Init(30000);
-   TTriangles1.Init(20000);
+}
+//=======================================================================
+//function : MakeSampling
+//purpose  :
+//=======================================================================
+void IntPolyh_MaillageAffinage::MakeSampling(const Standard_Integer SurfID,
+                                             TColStd_Array1OfReal& theUPars,
+                                             TColStd_Array1OfReal& theVPars)
+{
+  if (SurfID == 1)
+    IntPolyh_Tools::MakeSampling(MaSurface1, NbSamplesU1, NbSamplesV1, myEnlargeZone, theUPars, theVPars);
+  else
+    IntPolyh_Tools::MakeSampling(MaSurface2, NbSamplesU2, NbSamplesV2, myEnlargeZone, theUPars, theVPars);
+}
 
-   TPoints2.Init(10000);
-   TEdges2.Init(30000);
-   TTriangles2.Init(20000);
-   
-   TStartPoints.Init(10000);
- }
 //=======================================================================
 //function : FillArrayOfPnt
 //purpose  : Compute points on one surface and fill an array of points
@@ -218,74 +295,11 @@ IntPolyh_MaillageAffinage::IntPolyh_MaillageAffinage
 void IntPolyh_MaillageAffinage::FillArrayOfPnt
   (const Standard_Integer SurfID)
 {
-  Standard_Integer NbSamplesU, NbSamplesV;
-  Standard_Real u0, u1, v0, v1;
-  Handle(Adaptor3d_HSurface) MaSurface;
-  //
-  MaSurface=(SurfID==1)? MaSurface1:MaSurface2;
-  IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
-  NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
-  NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
-  //
-  u0 = (MaSurface)->FirstUParameter();  
-  u1 = (MaSurface)->LastUParameter();
-  v0 = (MaSurface)->FirstVParameter();  
-  v1 = (MaSurface)->LastVParameter();  
-
-  if(myEnlargeZone) {
-    if(MaSurface->GetType() == GeomAbs_BSplineSurface ||
-       MaSurface->GetType() == GeomAbs_BezierSurface) {
-      if((!MaSurface->IsUClosed() && !MaSurface->IsUPeriodic()) &&
-        (Abs(u0) < 1.e+100 && Abs(u1) < 1.e+100) ) {
-       Standard_Real delta_u = Abs(u1 - u0) / 100.;
-       u0 -= delta_u;
-       u1 += delta_u;
-      }
-      if((!MaSurface->IsVClosed() && !MaSurface->IsVPeriodic()) &&
-        (Abs(v0) < 1.e+100 && Abs(v1) < 1.e+100) ) {
-       Standard_Real delta_v = Abs(v1 - v0) / 100.;
-       v0 -= delta_v;
-       v1 += delta_v;
-      }
-    }
-  }
-  //
-  Standard_Integer iCnt, BoucleU, BoucleV;
-  Standard_Real itU, itV, U, V, Tol;
-  Bnd_Box *PtrBox;
-  gp_Pnt PtXYZ;
-  //
-  iCnt=0;
-  itU=(u1-u0)/Standard_Real(NbSamplesU-1);
-  itV=(v1-v0)/Standard_Real(NbSamplesV-1);
-  PtrBox = (SurfID==1) ? (&MyBox1) : (&MyBox2);
-
-  for(BoucleU=0; BoucleU<NbSamplesU; BoucleU++){
-     U = (BoucleU == (NbSamplesU - 1)) ? u1 : u0+BoucleU*itU;
-    for(BoucleV=0; BoucleV<NbSamplesV; BoucleV++){
-      V = (BoucleV == (NbSamplesV - 1)) ? v1 : v0+BoucleV*itV;
-      PtXYZ = (MaSurface)->Value(U,V);
-      IntPolyh_Point& aIPnt=TPoints[iCnt];
-      aIPnt.Set(PtXYZ.X(), PtXYZ.Y(), PtXYZ.Z(), U, V);
-      iCnt++;
-      PtrBox->Add(PtXYZ);
-    }
-  }
-  TPoints.SetNbPoints(iCnt);
-  //
-  IntCurveSurface_ThePolyhedronOfHInter polyhedron(MaSurface,
-                                                  NbSamplesU,
-                                                  NbSamplesV,
-                                                  u0,v0,
-                                                  u1,v1);
-  Tol=polyhedron.DeflectionOverEstimation();
-  Tol*=1.2;
-  //
-  Standard_Real a1,a2,a3,b1,b2,b3;
-  //
-  PtrBox->Get(a1,a2,a3,b1,b2,b3);
-  PtrBox->Update(a1-Tol,a2-Tol,a3-Tol,b1+Tol,b2+Tol,b3+Tol);
-  PtrBox->Enlarge(MyTolerance);
+  // Make sampling
+  TColStd_Array1OfReal aUpars, aVpars;
+  MakeSampling(SurfID, aUpars, aVpars);
+  // Fill array of points
+  FillArrayOfPnt(SurfID, aUpars, aVpars);
 }
 //=======================================================================
 //function : FillArrayOfPnt
@@ -296,86 +310,11 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
   (const Standard_Integer SurfID,
    const Standard_Boolean isShiftFwd)
 {
-
-  Handle(Adaptor3d_HSurface) MaSurface=(SurfID==1)? MaSurface1:MaSurface2;
-  IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
-  Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
-  Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
-
-  Standard_Real u0 = (MaSurface)->FirstUParameter();  
-  Standard_Real u1 = (MaSurface)->LastUParameter();
-  Standard_Real v0 = (MaSurface)->FirstVParameter();  
-  Standard_Real v1 = (MaSurface)->LastVParameter();  
-
-  if(myEnlargeZone) {
-    if(MaSurface->GetType() == GeomAbs_BSplineSurface ||
-       MaSurface->GetType() == GeomAbs_BezierSurface) {
-      if((!MaSurface->IsUClosed() && !MaSurface->IsUPeriodic()) &&
-        (Abs(u0) < 1.e+100 && Abs(u1) < 1.e+100) ) {
-       Standard_Real delta_u = Abs(u1 - u0) / 100.;
-       u0 -= delta_u;
-       u1 += delta_u;
-      }
-      if((!MaSurface->IsVClosed() && !MaSurface->IsVPeriodic()) &&
-        (Abs(v0) < 1.e+100 && Abs(v1) < 1.e+100) ) {
-       Standard_Real delta_v = Abs(v1 - v0) / 100.;
-       v0 -= delta_v;
-       v1 += delta_v;
-      }
-    }
-  }
-  
-  IntCurveSurface_ThePolyhedronOfHInter polyhedron(MaSurface,
-                                                  NbSamplesU,
-                                                  NbSamplesV,
-                                                  u0,v0,
-                                                  u1,v1);
-  Standard_Real Tol=polyhedron.DeflectionOverEstimation();
-
-  Standard_Integer CpteurTabPnt=0;
-  Standard_Real itU=(u1-u0)/Standard_Real(NbSamplesU-1);
-  Standard_Real itV=(v1-v0)/Standard_Real(NbSamplesV-1);
-
-  Bnd_Box *PtrBox = (SurfID==1) ? (&MyBox1) : (&MyBox2);
-  Standard_Real resol = gp::Resolution();
-
-  for(Standard_Integer BoucleU=0; BoucleU<NbSamplesU; BoucleU++){
-    Standard_Real U = (BoucleU == (NbSamplesU - 1)) ? u1 : u0+BoucleU*itU;
-    for(Standard_Integer BoucleV=0; BoucleV<NbSamplesV; BoucleV++){
-      Standard_Real V = (BoucleV == (NbSamplesV - 1)) ? v1 : v0+BoucleV*itV;
-
-      gp_Pnt PtXYZ;
-      gp_Vec aDU;
-      gp_Vec aDV;
-      gp_Vec aNorm;
-
-      MaSurface->D1(U, V, PtXYZ, aDU, aDV);
-      
-      aNorm = aDU.Crossed(aDV);
-      Standard_Real aMag = aNorm.Magnitude();
-      if (aMag > resol) {
-        aNorm /= aMag;
-        aNorm.Multiply(Tol*1.5);
-
-        if (isShiftFwd)
-          PtXYZ.Translate(aNorm);
-        else
-          PtXYZ.Translate(aNorm.Reversed());
-      }
-
-      (TPoints[CpteurTabPnt]).Set(PtXYZ.X(), PtXYZ.Y(), PtXYZ.Z(), U, V);
-      CpteurTabPnt++;
-      PtrBox->Add(PtXYZ);
-    }
-  }
-  TPoints.SetNbPoints(CpteurTabPnt);
-
-  Tol*=1.2;
-
-  Standard_Real a1,a2,a3,b1,b2,b3;
-  PtrBox->Get(a1,a2,a3,b1,b2,b3);
-  PtrBox->Update(a1-Tol,a2-Tol,a3-Tol,b1+Tol,b2+Tol,b3+Tol);
-  PtrBox->Enlarge(MyTolerance);
+  // Make sampling
+  TColStd_Array1OfReal aUpars, aVpars;
+  MakeSampling(SurfID, aUpars, aVpars);
+  // Fill array of points
+  FillArrayOfPnt(SurfID, isShiftFwd, aUpars, aVpars);
 }
 //=======================================================================
 //function : FillArrayOfPnt
@@ -384,7 +323,8 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
 void IntPolyh_MaillageAffinage::FillArrayOfPnt
   (const Standard_Integer SurfID, 
    const TColStd_Array1OfReal& Upars,
-   const TColStd_Array1OfReal& Vpars)
+   const TColStd_Array1OfReal& Vpars,
+   const Standard_Real *theDeflTol)
 {
   Standard_Boolean bDegI, bDeg;
   Standard_Integer aNbU, aNbV, iCnt, i, j;
@@ -398,7 +338,6 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
   Handle(Adaptor3d_HSurface)& aS=(SurfID==1)? MaSurface1:MaSurface2;
   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
   //
-  //modified by NIZNHY-PKV Fri Jan 20 09:48:57 2012f
   aJD1=0;
   aJD2=0;
   aID1=0;
@@ -407,13 +346,11 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
   if (!(aJD1 || aJD2)) {
     DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
   }
-  //modified by NIZNHY-PKV Fri Jan 20 09:49:00 2012t
   //
+  TPoints.Init(aNbU*aNbV);
   iCnt=0;
   for(i=1; i<=aNbU; ++i){
-    //modified by NIZNHY-PKV Fri Jan 20 13:59:15 2012f
     bDegI=(aID1==i || aID2==i);
-    //modified by NIZNHY-PKV Fri Jan 20 13:59:17 2012t
     aU=Upars(i);
     for(j=1; j<=aNbV; ++j){
       aV=Vpars(j);
@@ -422,23 +359,19 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
       IntPolyh_Point& aIP=TPoints[iCnt];
       aIP.Set(aX, aY, aZ, aU, aV);
       //
-      //modified by NIZNHY-PKV Fri Jan 20 13:59:06 2012f
       bDeg=bDegI || (aJD1==j || aJD2==j);
       if (bDeg) {
-       aIP.SetDegenerated(bDeg);
+        aIP.SetDegenerated(bDeg);
       }
-      //modified by NIZNHY-PKV Fri Jan 20 13:59:02 2012t
       ++iCnt;
       aBox.Add(aP);
     }
   }
   //
-  TPoints.SetNbPoints(iCnt);
+  TPoints.SetNbItems(iCnt);
   //
-  IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
-  //
-  aTol=polyhedron.DeflectionOverEstimation();
-  aTol*=1.2;
+  aTol = !theDeflTol ? IntPolyh_Tools::ComputeDeflection(aS, Upars, Vpars) : *theDeflTol;
+  aTol *= 1.2;
 
   Standard_Real a1,a2,a3,b1,b2,b3;
   //
@@ -449,97 +382,101 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
 
 //=======================================================================
 //function : FillArrayOfPnt
-//purpose  : Compute points on one surface and fill an array of points
-//           REMPLISSAGE DU TABLEAU DE POINTS
-//=======================================================================
-void IntPolyh_MaillageAffinage::FillArrayOfPnt
-  (const Standard_Integer SurfID,
-   const Standard_Boolean isShiftFwd,
-   const TColStd_Array1OfReal& Upars,
-   const TColStd_Array1OfReal& Vpars)
+//purpose  :
+//=======================================================================
+void IntPolyh_MaillageAffinage::FillArrayOfPnt(const Standard_Integer SurfID,
+                                               const Standard_Boolean isShiftFwd,
+                                               const IntPolyh_ArrayOfPointNormal& thePointsNorm,
+                                               const TColStd_Array1OfReal& theUPars,
+                                               const TColStd_Array1OfReal& theVPars,
+                                               const Standard_Real theDeflTol)
 {
-  Standard_Boolean bDegI, bDeg;
-  Standard_Integer aNbU, aNbV, iCnt, i, j;
-  Standard_Integer aID1, aID2, aJD1, aJD2;
-  Standard_Real Tol, resol, u0, v0, u1, v1, aU, aV, aMag;
-  Standard_Real aX, aY, aZ;
-  gp_Pnt aP;
-  gp_Vec aDU, aDV, aNorm;
-  //
-  aNbU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
-  aNbV=(SurfID==1)? NbSamplesV1:NbSamplesV2; 
+  Handle(Adaptor3d_HSurface) aS = (SurfID == 1) ? MaSurface1 : MaSurface2;
+  IntPolyh_ArrayOfPoints& TPoints = (SurfID == 1) ? TPoints1 : TPoints2;
+  Standard_Integer aNbU = (SurfID == 1) ? NbSamplesU1 : NbSamplesU2;
+  Standard_Integer aNbV = (SurfID == 1) ? NbSamplesV1 : NbSamplesV2;
   Bnd_Box& aBox = (SurfID==1) ? MyBox1 : MyBox2;
-  Handle(Adaptor3d_HSurface) aS=(SurfID==1)? MaSurface1:MaSurface2;
-  IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
-  //
-  resol = gp::Resolution();
-  u0 = Upars(1);
-  v0 = Vpars(1);
-  u1 = Upars(aNbU);
-  v1 = Vpars(aNbV);
-  //
-  IntCurveSurface_ThePolyhedronOfHInter polyhedron(aS, Upars, Vpars);
-  Tol=polyhedron.DeflectionOverEstimation();
-  //modified by NIZNHY-PKV Fri Jan 20 09:48:57 2012f
-  aJD1=0;
-  aJD2=0;
-  aID1=0;
-  aID2=0;
-  DegeneratedIndex(Vpars, aNbV, aS, 1, aJD1, aJD2);
-  if (!(aJD1 || aJD2)) {
-    DegeneratedIndex(Upars, aNbU, aS, 2, aID1, aID2);
-  }
-  //modified by NIZNHY-PKV Fri Jan 20 09:49:00 2012t
-  //
-  iCnt=0;
-  for(i=1; i<=aNbU; ++i){
-    //modified by NIZNHY-PKV Fri Jan 20 13:59:15 2012f
-    bDegI=(aID1==i || aID2==i);
-    //modified by NIZNHY-PKV Fri Jan 20 13:59:17 2012t
-    aU = Upars(i);
-    for(j=1; j<=aNbV; ++j){
-      aV = Vpars(j);
-      aS->D1(aU, aV, aP, aDU, aDV);
-      
-      aNorm = aDU.Crossed(aDV);
-      aMag = aNorm.Magnitude();
-      if (aMag > resol) {
-       aNorm /= aMag;
-       aNorm.Multiply(Tol*1.5);
-       //
-       if (isShiftFwd) {
-         aP.Translate(aNorm);
-       }
-       else{
-         aP.Translate(aNorm.Reversed());
-       }
-      }
-      
-      IntPolyh_Point& aIP=TPoints[iCnt];
+
+  Standard_Integer aJD1(0), aJD2(0), aID1(0), aID2(0);
+  DegeneratedIndex(theVPars, aNbV, aS, 1, aJD1, aJD2);
+  if (!(aJD1 || aJD2))
+    DegeneratedIndex(theUPars, aNbU, aS, 2, aID1, aID2);
+
+  Standard_Boolean bDegI, bDeg;
+  Standard_Integer iCnt(0), i, j;
+  Standard_Real aX, aY, aZ, aU, aV;
+
+  TPoints.Init(thePointsNorm.NbItems());
+
+  for (i = 1; i <= aNbU; ++i)
+  {
+    aU = theUPars(i);
+    bDegI = (aID1 == i || aID2 == i);
+    for (j = 1; j <= aNbV; ++j)
+    {
+      aV = theVPars(j);
+
+      const IntPolyh_PointNormal& aPN = thePointsNorm.Value(iCnt);
+      gp_Vec aNorm = aPN.Normal.Multiplied(1.5*theDeflTol);
+      if (!isShiftFwd)
+        aNorm.Reverse();
+      gp_Pnt aP = aPN.Point.Translated(aNorm);
+
+      IntPolyh_Point& aIP = TPoints[iCnt];
       aP.Coord(aX, aY, aZ);
       aIP.Set(aX, aY, aZ, aU, aV);
-      //
-      //modified by NIZNHY-PKV Fri Jan 20 13:59:06 2012f
-      bDeg=bDegI || (aJD1==j || aJD2==j);
-      if (bDeg) {
-       aIP.SetDegenerated(bDeg);
-      }
-      //modified by NIZNHY-PKV Fri Jan 20 13:59:02 2012t
+      bDeg = bDegI || (aJD1 == j || aJD2 == j);
+      if (bDeg)
+        aIP.SetDegenerated(bDeg);
+
       ++iCnt;
       aBox.Add(aP);
     }
   }
-  //
-  TPoints.SetNbPoints(iCnt);
-  //
-  Tol*=1.2;
-  //
+
+  TPoints.SetNbItems(iCnt);
+
+  // Update box
+  Standard_Real Tol = theDeflTol*1.2;
   Standard_Real a1,a2,a3,b1,b2,b3;
-  //
   aBox.Get(a1,a2,a3,b1,b2,b3);
   aBox.Update(a1-Tol,a2-Tol,a3-Tol,b1+Tol,b2+Tol,b3+Tol);
   aBox.Enlarge(MyTolerance);
 }
+
+//=======================================================================
+//function : FillArrayOfPnt
+//purpose  : Compute points on one surface and fill an array of points
+//=======================================================================
+void IntPolyh_MaillageAffinage::FillArrayOfPnt
+  (const Standard_Integer SurfID,
+   const Standard_Boolean isShiftFwd,
+   const TColStd_Array1OfReal& Upars,
+   const TColStd_Array1OfReal& Vpars,
+   const Standard_Real *theDeflTol)
+{
+  Handle(Adaptor3d_HSurface) aS = (SurfID == 1) ? MaSurface1 : MaSurface2;
+  // Compute the tolerance
+  Standard_Real aTol = theDeflTol != NULL ? * theDeflTol :
+    IntPolyh_Tools::ComputeDeflection(aS, Upars, Vpars);
+  // Fill array of point normal
+  IntPolyh_ArrayOfPointNormal aPoints;
+  IntPolyh_Tools::FillArrayOfPointNormal(aS, Upars, Vpars, aPoints);
+
+  // Fill array of points
+  FillArrayOfPnt(1, isShiftFwd, aPoints, Upars, Vpars, aTol);
+}
+
+//=======================================================================
+//function : CommonBox
+//purpose  : 
+//=======================================================================
+void IntPolyh_MaillageAffinage::CommonBox()
+{
+  Standard_Real XMin, YMin, ZMin, XMax, YMax, ZMax;
+  CommonBox(GetBox(1), GetBox(2), XMin, YMin, ZMin, XMax, YMax, ZMax);
+}
+
 //=======================================================================
 //function : CommonBox
 //purpose  : Compute the common box  witch is the intersection
@@ -549,13 +486,13 @@ void IntPolyh_MaillageAffinage::FillArrayOfPnt
 //           DETERMINATION OF THE COMMON BOX
 //=======================================================================
 void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
-                                          const Bnd_Box &,
-                                          Standard_Real &XMin,
-                                          Standard_Real &YMin,
-                                          Standard_Real &ZMin,
-                                          Standard_Real &XMax,
-                                          Standard_Real &YMax,
-                                          Standard_Real &ZMax) 
+                                           const Bnd_Box &,
+                                           Standard_Real &XMin,
+                                           Standard_Real &YMin,
+                                           Standard_Real &ZMin,
+                                           Standard_Real &XMax,
+                                           Standard_Real &YMax,
+                                           Standard_Real &ZMax) 
 {
   Standard_Real x10,y10,z10,x11,y11,z11;
   Standard_Real x20,y20,z20,x21,y21,z21;
@@ -580,8 +517,8 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     if(z20<=z10) ZMin=z10; else { if(z10<=z20) ZMin=z20;}
 
     if(((XMin==XMax)&&(!(YMin==YMax)&&!(ZMin==ZMax)))
-       ||((YMin==YMax)&&(!(XMin==XMax)&&!(ZMin==ZMax)))//ou exclusif ??
-       ||((ZMin==ZMax)&&(!(XMin==XMax)&&!(YMin==YMax)))) {
+        ||((YMin==YMax)&&(!(XMin==XMax)&&!(ZMin==ZMax)))//ou exclusif ??
+        ||((ZMin==ZMax)&&(!(XMin==XMax)&&!(YMin==YMax)))) {
     }
   }
   //
@@ -612,7 +549,7 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
   ZMin-=Z; ZMax+=Z;
 
   //Marking of points included in the common
-  const Standard_Integer FinTP1 = TPoints1.NbPoints();
+  const Standard_Integer FinTP1 = TPoints1.NbItems();
 //  for(Standard_Integer i=0; i<FinTP1; i++) {
   Standard_Integer i ;
   for( i=0; i<FinTP1; i++) {
@@ -623,9 +560,9 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     }   
     else { 
       if(Pt1.X()>XMax) {
-       r=2;
+        r=2;
       } else {
-       r=0; 
+        r=0; 
       } 
     }
     if(Pt1.Y()<YMin) { 
@@ -633,20 +570,20 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     }  
     else { 
       if(Pt1.Y()>YMax) {
-       r|=8; 
+        r|=8; 
       }
     }
     if(Pt1.Z()<ZMin) {
       r|=16; 
     } else {
       if(Pt1.Z()>ZMax) {
-       r|=32;
+        r|=32;
       }
     }
     Pt1.SetPartOfCommon(r);
   }
 
-  const Standard_Integer FinTP2 = TPoints2.NbPoints();
+  const Standard_Integer FinTP2 = TPoints2.NbItems();
   for(Standard_Integer ii=0; ii<FinTP2; ii++) {
     IntPolyh_Point & Pt2 = TPoints2[ii];
     Standard_Integer rr;
@@ -655,9 +592,9 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     }   
     else { 
       if(Pt2.X()>XMax) {
-       rr=2;
+        rr=2;
       } else {
-       rr=0;
+        rr=0;
       } 
     }
     if(Pt2.Y()<YMin) {
@@ -665,7 +602,7 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     }  
     else { 
       if(Pt2.Y()>YMax) {
-       rr|=8; 
+        rr|=8; 
       }
     }
     if(Pt2.Z()<ZMin) {
@@ -673,7 +610,7 @@ void IntPolyh_MaillageAffinage::CommonBox (const Bnd_Box &,
     }
     else {
       if(Pt2.Z()>ZMax) {
-       rr|=32;
+        rr|=32;
       }
     }
     Pt2.SetPartOfCommon(rr);
@@ -689,28 +626,37 @@ void IntPolyh_MaillageAffinage::FillArrayOfEdges
 {
 
   IntPolyh_ArrayOfEdges &TEdges=(SurfID==1)? TEdges1:TEdges2;
+  IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
 
+  //NbEdges = 3 + 3*(NbSamplesV-2) + 3*(NbSamplesU-2) + 
+  //        + 3*(NbSamplesU-2)*(NbSamplesV-2) + (NbSamplesV-1) + (NbSamplesU-1);
+  //NbSamplesU and NbSamples cannot be less than 2, so
+  Standard_Integer NbEdges = 3*NbSamplesU*NbSamplesV - 2*(NbSamplesU+NbSamplesV) + 1;
+  TEdges.Init(NbEdges);
+
   Standard_Integer CpteurTabEdges=0;
 
   //maillage u0 v0
   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
   TEdges[CpteurTabEdges].SetSecondPoint(1);             // U V+1
-  //  TEdges[CpteurTabEdges].SetFirstTriangle(-1);
   TEdges[CpteurTabEdges].SetSecondTriangle(0);
+  TTriangles[0].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
   CpteurTabEdges++;
-  
+
   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV);    // U+1 V
-  TEdges[CpteurTabEdges].SetFirstTriangle(0);
-  TEdges[CpteurTabEdges].SetSecondTriangle(1);
+  TEdges[CpteurTabEdges].SetFirstTriangle(1);
+  TTriangles[1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
   CpteurTabEdges++;
   
   TEdges[CpteurTabEdges].SetFirstPoint(0);                // U V
   TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV+1);  // U+1 V+1
-  TEdges[CpteurTabEdges].SetFirstTriangle(1);
-  //  TEdges[CpteurTabEdges].SetSecondTriangle(-1);
+  TEdges[CpteurTabEdges].SetFirstTriangle(0);
+  TTriangles[0].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
+  TEdges[CpteurTabEdges].SetSecondTriangle(1);
+  TTriangles[1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
   CpteurTabEdges++;
   
   //maillage surU=u0
@@ -721,41 +667,50 @@ void IntPolyh_MaillageAffinage::FillArrayOfEdges
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
     //    TEdges[CpteurTabEdges].SetFirstTriangle(-1);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2);
+    TTriangles[BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1
     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2);
+    TTriangles[BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2+1);
+    TTriangles[BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*2+1);
+    TTriangles[BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2-2);
+    TTriangles[BoucleMeshV*2-2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     PntInit++;
   }  
-       
+        
   //maillage sur V=v0
   PntInit=NbSamplesV;
   for(BoucleMeshV=1; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);    // U V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1); // U V+1
     TEdges[CpteurTabEdges].SetFirstTriangle((BoucleMeshV-1)*(NbSamplesV-1)*2+1);
+    TTriangles[(BoucleMeshV-1)*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2);
+    TTriangles[BoucleMeshV*(NbSamplesV-1)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit); // U V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);     // U+1 V+1
     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2);
+    TTriangles[BoucleMeshV*(NbSamplesV-1)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
+    TTriangles[BoucleMeshV*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);  // U V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);    // U+1 V
     TEdges[CpteurTabEdges].SetFirstTriangle(BoucleMeshV*(NbSamplesV-1)*2+1);
-    //    TEdges[CpteurTabEdges].SetSecondTriangle(-1);
+    TTriangles[BoucleMeshV*(NbSamplesV-1)*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     PntInit+=NbSamplesV;
   }
@@ -763,37 +718,43 @@ void IntPolyh_MaillageAffinage::FillArrayOfEdges
   PntInit=NbSamplesV+1;
   //To provide recursion I associate a point with three edges  
   for(Standard_Integer BoucleMeshU=1; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
-    for(Standard_Integer BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
+    for(BoucleMeshV=1; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);             // U V+1
       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*(BoucleMeshU-1)+BoucleMeshV*2+1);
+      TTriangles[(NbSamplesV-1)*2*(BoucleMeshU-1)+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
+      TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
       CpteurTabEdges++;
       
       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV+1);    // U+1 V+1 
       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2);
+      TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
+      TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
       CpteurTabEdges++;
       
       TEdges[CpteurTabEdges].SetFirstPoint(PntInit);                // U V
       TEdges[CpteurTabEdges].SetSecondPoint(PntInit+NbSamplesV);  // U+1 V
       TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1);
+      TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
       TEdges[CpteurTabEdges].SetSecondTriangle((NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2-2);
-      CpteurTabEdges++;            
+      TTriangles[(NbSamplesV-1)*2*BoucleMeshU+BoucleMeshV*2-2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
+      CpteurTabEdges++;            
       PntInit++;//Pass to the next point
     }
     PntInit++;//Pass the last point of the column
     PntInit++;//Pass the first point of the next column
   }
-       
+        
   //close mesh on U=u1
   PntInit=(NbSamplesU-1)*NbSamplesV; //point U=u1 V=0
   for(BoucleMeshV=0; BoucleMeshV<NbSamplesV-1; BoucleMeshV++){
     TEdges[CpteurTabEdges].SetFirstPoint(PntInit);           //U=u1 V
     TEdges[CpteurTabEdges].SetSecondPoint(PntInit+1);        //U=u1 V+1
     TEdges[CpteurTabEdges].SetFirstTriangle((NbSamplesU-2)*(NbSamplesV-1)*2+BoucleMeshV*2+1);
-    //    TEdges[CpteurTabEdges].SetSecondTriangle(-1);
+    TTriangles[(NbSamplesU-2)*(NbSamplesV-1)*2+BoucleMeshV*2+1].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
     PntInit++;
   }
@@ -802,12 +763,11 @@ void IntPolyh_MaillageAffinage::FillArrayOfEdges
   for(BoucleMeshV=0; BoucleMeshV<NbSamplesU-1;BoucleMeshV++){      
     TEdges[CpteurTabEdges].SetFirstPoint(NbSamplesV-1+BoucleMeshV*NbSamplesV);       // U V=v1
     TEdges[CpteurTabEdges].SetSecondPoint(NbSamplesV-1+(BoucleMeshV+1)*NbSamplesV);  //U+1 V=v1
-    //    TEdges[CpteurTabEdges].SetFirstTriangle(-1);
     TEdges[CpteurTabEdges].SetSecondTriangle(BoucleMeshV*2*(NbSamplesV-1)+(NbSamplesV-2)*2);
+    TTriangles[BoucleMeshV*2*(NbSamplesV-1)+(NbSamplesV-2)*2].SetEdgeAndOrientation(TEdges[CpteurTabEdges], CpteurTabEdges);
     CpteurTabEdges++;
   }
-  TEdges.SetNbEdges(CpteurTabEdges);
-
+  TEdges.SetNbItems(CpteurTabEdges);
 }
 
 //=======================================================================
@@ -828,7 +788,7 @@ void IntPolyh_MaillageAffinage::FillArrayOfTriangles
   Standard_Integer NbSamplesU=(SurfID==1)? NbSamplesU1:NbSamplesU2;
   Standard_Integer NbSamplesV=(SurfID==1)? NbSamplesV1:NbSamplesV2;
 
-  
+  TTriangles.Init(2*(NbSamplesU-1)*(NbSamplesV-1));
   //To provide recursion, I associate a point with two triangles  
   for(Standard_Integer BoucleMeshU=0; BoucleMeshU<NbSamplesU-1; BoucleMeshU++){
     for(Standard_Integer BoucleMeshV=0; BoucleMeshV<NbSamplesV-1;BoucleMeshV++){
@@ -840,10 +800,10 @@ void IntPolyh_MaillageAffinage::FillArrayOfTriangles
 
       // IF ITS EDGE CONTACTS WITH THE COMMON BOX IP REMAINS = A 1
       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+1].PartOfCommon()) )
-       &&( (TPoints[PntInit+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()))
-       &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
-       //IF NOT IP=0
-       TTriangles[CpteurTabTriangles].SetIndiceIntersectionPossible(0);
+         &&( (TPoints[PntInit+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()))
+         &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
+        //IF NOT IP=0
+        TTriangles[CpteurTabTriangles].SetIntersectionPossible(Standard_False);
 
       CpteurTabTriangles++;
 
@@ -854,9 +814,9 @@ void IntPolyh_MaillageAffinage::FillArrayOfTriangles
 
 
       if( ( (TPoints[PntInit].PartOfCommon()) & (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) )
-       &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV].PartOfCommon()))
-       &&( (TPoints[PntInit+NbSamplesV].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
-       TTriangles[CpteurTabTriangles].SetIndiceIntersectionPossible(0);
+         &&( (TPoints[PntInit+NbSamplesV+1].PartOfCommon()) & (TPoints[PntInit+NbSamplesV].PartOfCommon()))
+         &&( (TPoints[PntInit+NbSamplesV].PartOfCommon()) & (TPoints[PntInit].PartOfCommon())) ) 
+        TTriangles[CpteurTabTriangles].SetIntersectionPossible(Standard_False);
 
 
       CpteurTabTriangles++;
@@ -865,54 +825,27 @@ void IntPolyh_MaillageAffinage::FillArrayOfTriangles
     }
     PntInit++;//Pass the last point of the column
   }
-  TTriangles.SetNbTriangles(CpteurTabTriangles);
-  const Standard_Integer FinTT = TTriangles.NbTriangles();
+  TTriangles.SetNbItems(CpteurTabTriangles);
+  const Standard_Integer FinTT = TTriangles.NbItems();
   if (FinTT==0) {
   }
 }
 //=======================================================================
-//function : LinkEdges2Triangles
-//purpose  : fill the  edge fields in  Triangle object  for the
-//           two array of triangles.
-//=======================================================================
-void IntPolyh_MaillageAffinage::LinkEdges2Triangles() 
-{
-  const Standard_Integer FinTT1 = TTriangles1.NbTriangles();
-  const Standard_Integer FinTT2 = TTriangles2.NbTriangles();
-
-  for(Standard_Integer uiui1=0; uiui1<FinTT1; uiui1++) {
-    IntPolyh_Triangle & MyTriangle1=TTriangles1[uiui1];
-    if ( (MyTriangle1.FirstEdge()) == -1 ) {
-      MyTriangle1.SetEdgeandOrientation(1,TEdges1);
-      MyTriangle1.SetEdgeandOrientation(2,TEdges1);
-      MyTriangle1.SetEdgeandOrientation(3,TEdges1);
-    }
-  }
-  for(Standard_Integer uiui2=0; uiui2<FinTT2; uiui2++) {
-    IntPolyh_Triangle & MyTriangle2=TTriangles2[uiui2];
-    if ( (MyTriangle2.FirstEdge()) == -1 ) {
-      MyTriangle2.SetEdgeandOrientation(1,TEdges2);
-      MyTriangle2.SetEdgeandOrientation(2,TEdges2);
-      MyTriangle2.SetEdgeandOrientation(3,TEdges2);
-    }
-  }
-}
-//=======================================================================
 //function : CommonPartRefinement
 //purpose  : Refine systematicaly all marked triangles of both surfaces
 //           REFINING OF THE COMMON
 //=======================================================================
 void IntPolyh_MaillageAffinage::CommonPartRefinement() 
 {
-  Standard_Integer FinInit1 = TTriangles1.NbTriangles();
+  Standard_Integer FinInit1 = TTriangles1.NbItems();
   for(Standard_Integer i=0; i<FinInit1; i++) {
-    if(TTriangles1[i].IndiceIntersectionPossible()!=0) 
+    if(TTriangles1[i].IsIntersectionPossible())
       TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
   }
 
-  Standard_Integer FinInit2=TTriangles2.NbTriangles();
+  Standard_Integer FinInit2=TTriangles2.NbItems();
   for(Standard_Integer ii=0; ii<FinInit2; ii++) {
-    if(TTriangles2[ii].IndiceIntersectionPossible()!=0) 
+    if(TTriangles2[ii].IsIntersectionPossible())
       TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
   }
 
@@ -924,18 +857,18 @@ void IntPolyh_MaillageAffinage::CommonPartRefinement()
 void IntPolyh_MaillageAffinage::LocalSurfaceRefinement(const Standard_Integer SurfID) {
 //refine locally, but systematically the chosen surface
   if (SurfID==1) {
-    const Standard_Integer FinInit1 = TTriangles1.NbTriangles();
+    const Standard_Integer FinInit1 = TTriangles1.NbItems();
     for(Standard_Integer i=0; i<FinInit1; i++) {
-      if(TTriangles1[i].IndiceIntersectionPossible()!=0)
-       TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
+      if(TTriangles1[i].IsIntersectionPossible())
+        TTriangles1[i].MiddleRefinement(i,MaSurface1,TPoints1,TTriangles1,TEdges1);
     }
   }
   //
   if (SurfID==2) {
-    const Standard_Integer FinInit2 = TTriangles2.NbTriangles();
+    const Standard_Integer FinInit2 = TTriangles2.NbItems();
     for(Standard_Integer ii=0; ii<FinInit2; ii++) {
-      if(TTriangles2[ii].IndiceIntersectionPossible()!=0
-       TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
+      if(TTriangles2[ii].IsIntersectionPossible()
+        TTriangles2[ii].MiddleRefinement(ii,MaSurface2,TPoints2,TTriangles2,TEdges2); 
     }
   }
 }
@@ -951,391 +884,240 @@ void IntPolyh_MaillageAffinage::LocalSurfaceRefinement(const Standard_Integer Su
 void IntPolyh_MaillageAffinage::ComputeDeflections
   (const Standard_Integer SurfID)
 {
-  Handle(Adaptor3d_HSurface) MaSurface=(SurfID==1)? MaSurface1:MaSurface2;
+  Handle(Adaptor3d_HSurface) aSurface=(SurfID==1)? MaSurface1:MaSurface2;
   IntPolyh_ArrayOfPoints &TPoints=(SurfID==1)? TPoints1:TPoints2;
   IntPolyh_ArrayOfTriangles &TTriangles=(SurfID==1)? TTriangles1:TTriangles2;
   Standard_Real &FlecheMin=(SurfID==1)? FlecheMin1:FlecheMin2;
-  Standard_Real &FlecheMoy=(SurfID==1)? FlecheMoy1:FlecheMoy2;
   Standard_Real &FlecheMax=(SurfID==1)? FlecheMax1:FlecheMax2;
 
-  Standard_Integer CpteurTabFleche=0;
   FlecheMax=-RealLast();
   FlecheMin=RealLast();
-  FlecheMoy=0.0;
-  const Standard_Integer FinTT = TTriangles.NbTriangles();
+  const Standard_Integer FinTT = TTriangles.NbItems();
   
-  for(CpteurTabFleche=0; CpteurTabFleche<FinTT; CpteurTabFleche++) {
-    IntPolyh_Triangle &Triangle = TTriangles[CpteurTabFleche];
-    if ( Triangle.GetFleche() < 0) { //pas normal
+  for(Standard_Integer i = 0; i < FinTT; i++) {
+    IntPolyh_Triangle& aTriangle = TTriangles[i];
+    Standard_Real Fleche = aTriangle.ComputeDeflection(aSurface, TPoints);
+    if (Fleche > FlecheMax)
+      FlecheMax = Fleche;
+    if (Fleche < FlecheMin)
+      FlecheMin = Fleche;
+  }
+}
 
+//=======================================================================
+//function : TrianglesDeflectionsRefinement
+//purpose  : Refinement of the triangles depending on the deflection
+//=======================================================================
+static
+  void TrianglesDeflectionsRefinement(const Handle(Adaptor3d_HSurface)& theS1,
+                                      IntPolyh_ArrayOfTriangles& theTriangles1,
+                                      IntPolyh_ArrayOfEdges& theEdges1,
+                                      IntPolyh_ArrayOfPoints& thePoints1,
+                                      const Standard_Real theFlecheCritique1,
+                                      const Handle(Adaptor3d_HSurface)& theS2,
+                                      IntPolyh_ArrayOfTriangles& theTriangles2,
+                                      IntPolyh_ArrayOfEdges& theEdges2,
+                                      IntPolyh_ArrayOfPoints& thePoints2,
+                                      const Standard_Real theFlecheCritique2)
+{
+  // Find the intersecting triangles
+  IntPolyh_IndexedDataMapOfIntegerArrayOfInteger aDMILI;
+  GetInterferingTriangles(theTriangles1, thePoints1, theTriangles2, thePoints2, aDMILI);
+  //
+  // Interfering triangles of second surface
+  TColStd_MapOfInteger aMIntS2;
+  // Since the number of the triangles may change during refinement,
+  // save the number of triangles before refinement
+  Standard_Integer FinTT1 = theTriangles1.NbItems();
+  Standard_Integer FinTT2 = theTriangles2.NbItems();
+  //
+  // Analyze interfering triangles
+  for (Standard_Integer i_S1 = 0; i_S1 < FinTT1; i_S1++) {
+    IntPolyh_Triangle& aTriangle1 = theTriangles1[i_S1];
+    if (!aTriangle1.IsIntersectionPossible()) {
+      continue;
+    }
+    //
+    const IntPolyh_ArrayOfInteger *pLI = aDMILI.Seek(i_S1);
+    if (!pLI || !pLI->Length()) {
+      // Mark non-interfering triangles of S1 to avoid their repeated usage
+      aTriangle1.SetIntersectionPossible(Standard_False);
+      continue;
+    }
+    //
+    if (aTriangle1.Deflection() > theFlecheCritique1) {
+      aTriangle1.MiddleRefinement(i_S1, theS1, thePoints1, theTriangles1, theEdges1);
+    }
+    //
+    IntPolyh_ArrayOfInteger::Iterator Iter(*pLI);
+    for (; Iter.More(); Iter.Next()) {
+      Standard_Integer i_S2 = Iter.Value();
+      if (aMIntS2.Add(i_S2)) {
+        IntPolyh_Triangle & aTriangle2 = theTriangles2[i_S2];
+        if (aTriangle2.Deflection() > theFlecheCritique2) {
+          // Refinement of the larger triangles
+          aTriangle2.MiddleRefinement(i_S2, theS2, thePoints2, theTriangles2, theEdges2);
+        }
+      }
     }
-    else{
-      Triangle.TriangleDeflection(MaSurface, TPoints);
-      Standard_Real Fleche=Triangle.GetFleche();
-
-      if (Fleche > FlecheMax)
-       FlecheMax=Fleche;
-      if (Fleche < FlecheMin)
-       FlecheMin=Fleche;
+  }
+  //
+  // Mark non-interfering triangles of S2 to avoid their repeated usage
+  if (aMIntS2.Extent() < FinTT2) {
+    for (Standard_Integer i_S2 = 0; i_S2 < FinTT2; i_S2++) {
+      if (!aMIntS2.Contains(i_S2)) {
+        theTriangles2[i_S2].SetIntersectionPossible(Standard_False);
+      }
+    }
+  }
+}
+//=======================================================================
+//function : LargeTrianglesDeflectionsRefinement
+//purpose  : Refinement of the large triangles in case one surface is
+//           much smaller then the other.
+//=======================================================================
+static
+  void LargeTrianglesDeflectionsRefinement(const Handle(Adaptor3d_HSurface)& theS,
+                                           IntPolyh_ArrayOfTriangles& theTriangles,
+                                           IntPolyh_ArrayOfEdges& theEdges,
+                                           IntPolyh_ArrayOfPoints& thePoints,
+                                           const Bnd_Box& theOppositeBox)
+{
+  // Find all triangles of the bigger surface with bounding boxes
+  // overlapping the bounding box the other surface
+  TColStd_ListOfInteger aLIT;
+  Standard_Integer i, aNbT = theTriangles.NbItems();
+  for (i = 0; i < aNbT; ++i) {
+    IntPolyh_Triangle& aTriangle = theTriangles[i];
+    if (!aTriangle.IsIntersectionPossible() || aTriangle.IsDegenerated()) {
+      continue;
+    }
+    //
+    const Bnd_Box& aBox = aTriangle.BoundingBox(thePoints);
+    if (aBox.IsVoid()) {
+      continue;
+    }
+    //
+    if (aBox.IsOut(theOppositeBox)) {
+      aTriangle.SetIntersectionPossible(Standard_False);
+      continue;
+    }
+    //
+    aLIT.Append(i);
+  }
+  //
+  if (aLIT.IsEmpty()) {
+    return;
+  }
+  //
+  // One surface is very small relatively to the other.
+  // The criterion of refining for large surface depends on the size of
+  // the bounding box of the other - since the criterion should be minimized,
+  // the smallest side of the bounding box is taken
+  Standard_Real x0, y0, z0, x1, y1, z1;
+  theOppositeBox.Get(x0, y0, z0, x1, y1, z1);
+  Standard_Real dx = Abs(x1 - x0);
+  Standard_Real dy = Abs(y1 - y0);
+  Standard_Real diag = Abs(z1 - z0);
+  Standard_Real dd = (dx > dy) ? dy : dx;
+  if (diag > dd) diag = dd;
+
+  // calculation of the criterion of refining
+  Standard_Real CritereAffinage = 0.0;
+  Standard_Real DiagPonderation = 0.5;
+  CritereAffinage = diag*DiagPonderation;
+  //
+  // The deflection of the greater is compared to the size of the smaller
+  TColStd_ListIteratorOfListOfInteger Iter(aLIT);
+  for (; Iter.More(); Iter.Next()) {
+    i = Iter.Value();
+    IntPolyh_Triangle & aTriangle = theTriangles[i];
+    if (aTriangle.Deflection() > CritereAffinage) {
+      aTriangle.MultipleMiddleRefinement(CritereAffinage, theOppositeBox, i,
+                                         theS, thePoints, theTriangles, theEdges);
+    }
+    else {
+      aTriangle.MiddleRefinement(i, theS, thePoints, theTriangles, theEdges);
     }
   }
 }
 //=======================================================================
 //function : TrianglesDeflectionsRefinementBSB
-//purpose  : Refine  both  surfaces using  BoundSortBox  as --
-//           rejection.  The  criterions  used to refine a  --
-//           triangle are:  The deflection The  size of the --
-//           bounding boxes   (one surface may be   very small
-//           compared to the other)
+//purpose  : Refine both surfaces using UBTree as rejection.
+//           The criterion used to refine a triangle are:
+//           - The deflection;
+//           - The  size of the - bounding boxes (one surface may be
+//           very small compared to the other).
 //=======================================================================
-void IntPolyh_MaillageAffinage::TrianglesDeflectionsRefinementBSB() 
+void IntPolyh_MaillageAffinage::TrianglesDeflectionsRefinementBSB()
 {
-  const Standard_Integer FinTT1 = TTriangles1.NbTriangles();
-  const Standard_Integer FinTT2 = TTriangles2.NbTriangles();
-  
+  // To estimate a surface in general it can be interesting
+  // to calculate all deflections
   ComputeDeflections(1);
-  // To estimate a surface in general it can be interesting 
-  //to calculate all deflections
-  //-- Check deflection at output
-                                                                
+  // Check deflection at output
   Standard_Real FlecheCritique1;
-  if(FlecheMin1>FlecheMax1) {
-    return;    
+  if (FlecheMin1 > FlecheMax1) {
+    return;
   }
-  else {//fleche min + (flechemax-flechemin) * 80/100
-    FlecheCritique1 = FlecheMin1*0.2+FlecheMax1*0.8;
+  else {
+    // fleche min + (flechemax-flechemin) * 80/100
+    FlecheCritique1 = FlecheMin1*0.2 + FlecheMax1*0.8;
   }
-  
+
+  // Compute deflections for the second surface
   ComputeDeflections(2);
+
   //-- Check arrows at output
-  
   Standard_Real FlecheCritique2;
-  if(FlecheMin2>FlecheMax2) { 
-
+  if (FlecheMin2 > FlecheMax2) {
     return;
   }
-  else {//fleche min + (flechemax-flechemin) * 80/100
-    FlecheCritique2 = FlecheMin2*0.2+FlecheMax2*0.8;
-  }
-  
-  //Bounding boxes
-  Bnd_BoundSortBox BndBSB;
-  Standard_Real diag1,diag2;
-  Standard_Real x0,y0,z0,x1,y1,z1;
-  
-  //The greatest of two bounding boxes created in FillArrayOfPoints is found.
-  //Then this value is weighted depending on the discretization 
-  //(NbSamplesU and NbSamplesV)
-  MyBox1.Get(x0,y0,z0,x1,y1,z1);
-  x0-=x1; y0-=y1; z0-=z1;
-  diag1=x0*x0+y0*y0+z0*z0;
-  const Standard_Real NbSamplesUV1=Standard_Real(NbSamplesU1) * Standard_Real(NbSamplesV1);
-  diag1/=NbSamplesUV1;
-
-  MyBox2.Get(x0,y0,z0,x1,y1,z1);
-  x0-=x1; y0-=y1; z0-=z1;
-  diag2=x0*x0+y0*y0+z0*z0;
-  const Standard_Real NbSamplesUV2=Standard_Real(NbSamplesU2) * Standard_Real(NbSamplesV2);
-  diag2/=NbSamplesUV2;
-  
-  //-- The surface with the greatest bounding box is "discretized"
-  
-  //Standard_Integer NbInterTentees=0;
-  
-  if(diag1<diag2) { 
-    
-    if(FlecheCritique2<diag1) {//the corresponding sizes are not too disproportional
-
-      Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT2);
-        
-      for(Standard_Integer i=0; i<FinTT2; i++){
-       if (TTriangles2[i].IndiceIntersectionPossible()!=0) {
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles2[i];
-         const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()];
-         const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);//Box b, which contains triangle i of surface 2 is created./
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche()+MyTolerance);
-         HBnd->SetValue(i+1,b);//Box b is added in the array HBnd
-       }
-      }
-      
-      //Inititalization of the boundary, sorting of boxes
-      BndBSB.Initialize(HBnd);//contains boxes of 2
-      
-      Standard_Integer FinTT1Init=FinTT1;
-      for(Standard_Integer i_S1=0; i_S1<FinTT1Init; i_S1++) {
-       if(TTriangles1[i_S1].IndiceIntersectionPossible()!=0) {
-         //-- Loop on the boxes of mesh 1 
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles1[i_S1];
-         const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()]; 
-         const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche());
-         //-- List of boxes of 2, which touch this box (of 1)
-         const TColStd_ListOfInteger& ListeOf2 = BndBSB.Compare(b);
-         
-         if((ListeOf2.IsEmpty())==0) {
-           IntPolyh_Triangle &Triangle1 = TTriangles1[i_S1];
-           if(Triangle1.GetFleche()>FlecheCritique1)
-             Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
-                                        TTriangles1, TEdges1);
-           
-           for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf2); 
-                Iter.More(); 
-                Iter.Next()) {
-             Standard_Integer i_S2=Iter.Value()-1;
-             //if the box of s1 contacts with the boxes of s2 
-             //the arrow of the triangle is checked
-             IntPolyh_Triangle & Triangle2 = TTriangles2[i_S2];
-             if(Triangle2.IndiceIntersectionPossible()!=0)
-               if(Triangle2.GetFleche()>FlecheCritique2)
-                 Triangle2.MiddleRefinement( i_S2, MaSurface2, TPoints2,
-                                            TTriangles2, TEdges2);
-           }
-         }
-       }
-      }
-    }
-
-    //--------------------------------------------------------------------
-    //FlecheCritique2 > diag1
+  else {
+    //fleche min + (flechemax-flechemin) * 80/100
+    FlecheCritique2 = FlecheMin2*0.2 + FlecheMax2*0.8;
+  }
+
+  // The greatest of two bounding boxes created in FillArrayOfPoints is found.
+  // Then this value is weighted depending on the discretization 
+  // (NbSamplesU and NbSamplesV)
+  Standard_Real diag1, diag2;
+  Standard_Real x0, y0, z0, x1, y1, z1;
+
+  MyBox1.Get(x0, y0, z0, x1, y1, z1);
+  x0 -= x1; y0 -= y1; z0 -= z1;
+  diag1 = x0*x0 + y0*y0 + z0*z0;
+  const Standard_Real NbSamplesUV1 = Standard_Real(NbSamplesU1) * Standard_Real(NbSamplesV1);
+  diag1 /= NbSamplesUV1;
+
+  MyBox2.Get(x0, y0, z0, x1, y1, z1);
+  x0 -= x1; y0 -= y1; z0 -= z1;
+  diag2 = x0*x0 + y0*y0 + z0*z0;
+  const Standard_Real NbSamplesUV2 = Standard_Real(NbSamplesU2) * Standard_Real(NbSamplesV2);
+  diag2 /= NbSamplesUV2;
+
+  // The surface with the greatest bounding box is "discretized"
+  if (diag1 < diag2) {
+    // second is discretized
+    if (FlecheCritique2 < diag1) {
+      // The corresponding sizes are not too disproportional
+      TrianglesDeflectionsRefinement(MaSurface1, TTriangles1, TEdges1, TPoints1, FlecheCritique1,
+                                     MaSurface2, TTriangles2, TEdges2, TPoints2, FlecheCritique2);
+    }
     else {
-      //2 is discretized
-
-      Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT2);
-    
-      for(Standard_Integer i=0; i<FinTT2; i++){
-       if (TTriangles2[i].IndiceIntersectionPossible()!=0) {
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles2[i];
-         const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()];
-         const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);//Box b, which contains triangle i of surface 2 is created/
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche()+MyTolerance);
-         //-- BndBSB.Add(b,i+1);
-         HBnd->SetValue(i+1,b);//Box b is added in array HBnd
-       }
-      }
-      
-      //Inititalization of the ouput bounding box
-      BndBSB.Initialize(HBnd);//contains boxes of 2
-      
-
-      //The bounding box Be1 of surface1 is compared BSB of surface2
-      const TColStd_ListOfInteger& ListeOf2 = BndBSB.Compare(MyBox1);
-      
-      if((ListeOf2.IsEmpty())==0) {
-       //if the bounding box Be1 of s1 contacts with 
-       //the boxes of s2 the deflection of triangle of s2 is checked
-
-       // Be1 is very small in relation to Be2
-       //The criterion of refining for surface2 depends on the size of Be1
-       //As it is known that this criterion should be minimized, 
-       //the smallest side of the bounding box is taken
-       Standard_Real x0,x1,y0,y1,z0,z1;
-       MyBox1.Get(x0,y0,z0,x1,y1,z1);
-       Standard_Real dx=Abs(x1-x0);
-       Standard_Real dy=Abs(y1-y0);
-       Standard_Real diag=Abs(z1-z0);
-       Standard_Real dd=-1.0;
-       if (dx>dy)
-         dd=dy;
-       else
-         dd=dx;
-       if (diag>dd) diag=dd;
-       
-       //if Be1 contacts with the boxes of s2, the deflection 
-       //of the triangles of s2 is checked (greater)
-       //in relation to the size of Be1 (smaller)
-       for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf2); 
-            Iter.More(); 
-            Iter.Next()) {
-         Standard_Integer i_S2=Iter.Value()-1;
-         
-         IntPolyh_Triangle & Triangle2=TTriangles2[i_S2];
-         if(Triangle2.IndiceIntersectionPossible()) {
-           
-           //calculation of the criterion of refining
-           //The deflection of the greater is compared to the size of the smaller
-           Standard_Real CritereAffinage=0.0;
-           Standard_Real DiagPonderation=0.5;
-           CritereAffinage = diag*DiagPonderation;
-           if(Triangle2.GetFleche()>CritereAffinage)
-             Triangle2.MultipleMiddleRefinement2(CritereAffinage, MyBox1, i_S2,
-                                                 MaSurface2, TPoints2,
-                                                 TTriangles2,TEdges2);
-           
-           else Triangle2.MiddleRefinement(i_S2,MaSurface2,TPoints2,
-                                           TTriangles2, TEdges2);
-         }
-       }
-      }
+      // second surface is much larger then the first
+      LargeTrianglesDeflectionsRefinement(MaSurface2, TTriangles2, TEdges2, TPoints2, MyBox1);
     }
   }
-  
-  
-  else {     //-- The greater is discretised
-
-    if(FlecheCritique1<diag2) {//the respective sizes are not to much disproportional
-    
-      Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT1);
-      
-      for(Standard_Integer i=0; i<FinTT1; i++){
-       if(TTriangles1[i].IndiceIntersectionPossible()!=0) {
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles1[i];
-         const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()];
-         const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);//Box b, which contains triangle i of surface 2 is created.
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche()+MyTolerance);
-         HBnd->SetValue(i+1,b);//Boite b is added in the array HBnd
-       }
-      }
-      BndBSB.Initialize(HBnd);
-      
-      Standard_Integer FinTT2init=FinTT2;
-      for(Standard_Integer i_S2=0; i_S2<FinTT2init; i_S2++) {
-       if (TTriangles2[i_S2].IndiceIntersectionPossible()!=0) {      
-         //-- Loop on the boxes of mesh 2 
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles2[i_S2];
-         const IntPolyh_Point&    PA=TPoints2[T.FirstPoint()]; 
-         const IntPolyh_Point&    PB=TPoints2[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints2[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche()+MyTolerance);
-         //-- List of boxes of 1 touching this box (of 2)
-         const TColStd_ListOfInteger& ListeOf1 = BndBSB.Compare(b);
-         IntPolyh_Triangle & Triangle2=TTriangles2[i_S2];
-         if((ListeOf1.IsEmpty())==0) {
-           
-           if(Triangle2.GetFleche()>FlecheCritique2)
-             Triangle2.MiddleRefinement(i_S2,MaSurface2,TPoints2,
-                                        TTriangles2, TEdges2);
-           
-           for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf1); 
-                Iter.More(); 
-                Iter.Next()) {
-             Standard_Integer i_S1=Iter.Value()-1;
-             IntPolyh_Triangle & Triangle1=TTriangles1[i_S1];
-             if (Triangle1.IndiceIntersectionPossible())
-               if(Triangle1.GetFleche()>FlecheCritique1)
-                 Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
-                                            TTriangles1, TEdges1); 
-           }
-         }
-       }
-      }
-    }
-    //-----------------------------------------------------------------------------
-    else {// FlecheCritique1>diag2
-      // 1 is discretized
-
-      Handle(Bnd_HArray1OfBox) HBnd = new  Bnd_HArray1OfBox(1,FinTT1);
-    
-      for(Standard_Integer i=0; i<FinTT1; i++){
-       if (TTriangles1[i].IndiceIntersectionPossible()!=0) {
-         Bnd_Box b;
-         const IntPolyh_Triangle& T=TTriangles1[i];
-         const IntPolyh_Point&    PA=TPoints1[T.FirstPoint()];
-         const IntPolyh_Point&    PB=TPoints1[T.SecondPoint()]; 
-         const IntPolyh_Point&    PC=TPoints1[T.ThirdPoint()]; 
-         gp_Pnt pntA(PA.X(),PA.Y(),PA.Z());
-         gp_Pnt pntB(PB.X(),PB.Y(),PB.Z());
-         gp_Pnt pntC(PC.X(),PC.Y(),PC.Z());
-         b.Add(pntA);//Box b, which contains triangle i of surface 1 is created./
-         b.Add(pntB);
-         b.Add(pntC);
-         b.Enlarge(T.GetFleche()+MyTolerance);
-         HBnd->SetValue(i+1,b);//Box b is added in the array HBnd
-       }
-      }
-      
-      //Inititalisation of the boundary output box
-      BndBSB.Initialize(HBnd);//contains boxes of 1
-      
-      //Bounding box Be2 of surface2 is compared to BSB of surface1
-      const TColStd_ListOfInteger& ListeOf1 = BndBSB.Compare(MyBox2);
-      
-      if((ListeOf1.IsEmpty())==0) {
-       //if the bounding box Be2 of s2 contacts 
-       //with boxes of s1 the deflection of the triangle of s1 is checked
-       
-       // Be2 is very small compared to Be1
-       //The criterion of refining for surface1 depends on the size of Be2
-       //As this criterion should be minimized, 
-       //the smallest side of the bounding box is taken
-       Standard_Real x0,x1,y0,y1,z0,z1;
-       MyBox2.Get(x0,y0,z0,x1,y1,z1);
-       Standard_Real dx=Abs(x1-x0);
-       Standard_Real dy=Abs(y1-y0);
-       Standard_Real diag=Abs(z1-z0);
-       Standard_Real dd=-1.0;
-       if (dx>dy)
-         dd=dy;
-       else
-         dd=dx;
-       if (diag>dd) diag=dd;
-       
-       //if Be2 contacts with boxes of s1, the deflection of 
-       //triangles of s1 (greater) is checked
-       //comparatively to the size of Be2 (smaller).
-       for (TColStd_ListIteratorOfListOfInteger Iter(ListeOf1); 
-            Iter.More(); 
-            Iter.Next()) {
-         Standard_Integer i_S1=Iter.Value()-1;
-
-         IntPolyh_Triangle & Triangle1=TTriangles1[i_S1];
-         if(Triangle1.IndiceIntersectionPossible()) {
-
-           //calculation of the criterion of refining
-           //The deflection of the greater is compared 
-           //with the size of the smaller.
-           Standard_Real CritereAffinage=0.0;
-           Standard_Real DiagPonderation=0.5;
-           CritereAffinage = diag*DiagPonderation;;
-           if(Triangle1.GetFleche()>CritereAffinage)
-             Triangle1.MultipleMiddleRefinement2(CritereAffinage,MyBox2, i_S1,
-                                                 MaSurface1, TPoints1,
-                                                 TTriangles1, TEdges1); 
-           
-           else Triangle1.MiddleRefinement(i_S1,MaSurface1,TPoints1,
-                                           TTriangles1, TEdges1);
-           
-         }
-       }
-      }
+  else {
+    // first is discretized
+    if (FlecheCritique1 < diag2) {
+      // The corresponding sizes are not too disproportional
+      TrianglesDeflectionsRefinement(MaSurface2, TTriangles2, TEdges2, TPoints2, FlecheCritique2,
+                                     MaSurface1, TTriangles1, TEdges1, TPoints1, FlecheCritique1);
+    }
+    else {
+      // first surface is much larger then the second
+      LargeTrianglesDeflectionsRefinement(MaSurface1, TTriangles1, TEdges1, TPoints1, MyBox2);
     }
   }
 }
@@ -1344,8 +1126,8 @@ void IntPolyh_MaillageAffinage::TrianglesDeflectionsRefinementBSB()
 //purpose  : This function is used for the function project6
 //=======================================================================
 inline Standard_Real maxSR(const Standard_Real a,
-                          const Standard_Real b,
-                          const Standard_Real c)
+                           const Standard_Real b,
+                           const Standard_Real c)
 {
   Standard_Real t = a;
   if (b > t) t = b;
@@ -1357,8 +1139,8 @@ inline Standard_Real maxSR(const Standard_Real a,
 //purpose  : This function is used for the function project6
 //=======================================================================
 inline Standard_Real minSR(const Standard_Real a,
-                          const Standard_Real b,
-                          const Standard_Real c)
+                           const Standard_Real b,
+                           const Standard_Real c)
 {
   Standard_Real t = a;
   if (b < t) t = b;
@@ -1371,11 +1153,11 @@ inline Standard_Real minSR(const Standard_Real a,
 //=======================================================================
 Standard_Integer project6(const IntPolyh_Point &ax, 
                           const IntPolyh_Point &p1,
-                         const IntPolyh_Point &p2, 
-                         const IntPolyh_Point &p3, 
+                          const IntPolyh_Point &p2, 
+                          const IntPolyh_Point &p3, 
                           const IntPolyh_Point &q1,
-                         const IntPolyh_Point &q2, 
-                         const IntPolyh_Point &q3)
+                          const IntPolyh_Point &q2, 
+                          const IntPolyh_Point &q3)
 {
   Standard_Real P1 = ax.Dot(p1);
   Standard_Real P2 = ax.Dot(p2);
@@ -1413,7 +1195,15 @@ Standard_Integer IntPolyh_MaillageAffinage::TriContact
      The edges are (e1,e2,e3) and (f1,f2,f3).
      The normals are n1 and m1
      Outwards are (g1,g2,g3) and (h1,h2,h3).*/
+  
+  if(maxSR(P1.X(),P2.X(),P3.X())<minSR(Q1.X(),Q2.X(),Q3.X())) return(0);
+  if(maxSR(P1.Y(),P2.Y(),P3.Y())<minSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
+  if(maxSR(P1.Z(),P2.Z(),P3.Z())<minSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
 
+  if(minSR(P1.X(),P2.X(),P3.X())>maxSR(Q1.X(),Q2.X(),Q3.X())) return(0);
+  if(minSR(P1.Y(),P2.Y(),P3.Y())>maxSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
+  if(minSR(P1.Z(),P2.Z(),P3.Z())>maxSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
+    
   IntPolyh_Point p1, p2, p3;
   IntPolyh_Point q1, q2, q3;
   IntPolyh_Point e1, e2, e3;
@@ -1426,16 +1216,9 @@ Standard_Integer IntPolyh_MaillageAffinage::TriContact
   IntPolyh_Point ef11, ef12, ef13;
   IntPolyh_Point ef21, ef22, ef23;
   IntPolyh_Point ef31, ef32, ef33;
-  
+
   z.SetX(0.0);  z.SetY(0.0);  z.SetZ(0.0);
-  
-  if(maxSR(P1.X(),P2.X(),P3.X())<minSR(Q1.X(),Q2.X(),Q3.X())) return(0);
-  if(maxSR(P1.Y(),P2.Y(),P3.Y())<minSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
-  if(maxSR(P1.Z(),P2.Z(),P3.Z())<minSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
 
-  if(minSR(P1.X(),P2.X(),P3.X())>maxSR(Q1.X(),Q2.X(),Q3.X())) return(0);
-  if(minSR(P1.Y(),P2.Y(),P3.Y())>maxSR(Q1.Y(),Q2.Y(),Q3.Y())) return(0);
-  if(minSR(P1.Z(),P2.Z(),P3.Z())>maxSR(Q1.Z(),Q2.Z(),Q3.Z())) return(0);
 
   p1.SetX(P1.X() - P1.X());  p1.SetY(P1.Y() - P1.Y());  p1.SetZ(P1.Z() - P1.Z());
   p2.SetX(P2.X() - P1.X());  p2.SetY(P2.Y() - P1.Y());  p2.SetZ(P2.Z() - P1.Z());
@@ -1514,14 +1297,14 @@ Standard_Integer IntPolyh_MaillageAffinage::TriContact
 //           void TestNbPoints(const Standard_Integer TriSurfID,
 //=======================================================================
 void TestNbPoints(const Standard_Integer ,
-                 Standard_Integer &NbPoints,
-                 Standard_Integer &NbPointsTotal,
-                 const IntPolyh_StartPoint &Pt1,
-                 const IntPolyh_StartPoint &Pt2,
-                 IntPolyh_StartPoint &SP1,
-                 IntPolyh_StartPoint &SP2)
+                  Standard_Integer &NbPoints,
+                  Standard_Integer &NbPointsTotal,
+                  const IntPolyh_StartPoint &Pt1,
+                  const IntPolyh_StartPoint &Pt2,
+                  IntPolyh_StartPoint &SP1,
+                  IntPolyh_StartPoint &SP2)
 {
-  // already checked in TriangleEdgeContact2
+  // already checked in TriangleEdgeContact
   //  if( (NbPoints==2)&&(Pt1.CheckSameSP(Pt2)) ) NbPoints=1;
 
   if(NbPoints>2) {
@@ -1534,13 +1317,13 @@ void TestNbPoints(const Standard_Integer ,
     }
     else if ( (NbPoints==1)&&(NbPointsTotal==1) ) { 
       if(Pt1.CheckSameSP(SP1)!=1) {
-       SP2=Pt1;
-       NbPointsTotal=2;
+        SP2=Pt1;
+        NbPointsTotal=2;
       }
     }
     else if( (NbPoints==1)&&(NbPointsTotal==2) ) {
       if ( (SP1.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt1)) ) 
-       NbPointsTotal=2;
+        NbPointsTotal=2;
       else NbPointsTotal=3;
     }
     else if( (NbPoints==2)&&(NbPointsTotal==0) ) {
@@ -1550,20 +1333,20 @@ void TestNbPoints(const Standard_Integer ,
     }
     else if( (NbPoints==2)&&(NbPointsTotal==1) ) {//there is also Pt1 != Pt2 
       if(SP1.CheckSameSP(Pt1)) {
-       SP2=Pt2;
-       NbPointsTotal=2;
+        SP2=Pt2;
+        NbPointsTotal=2;
       }
       else if (SP1.CheckSameSP(Pt2)) {
-       SP2=Pt1;
-       NbPointsTotal=2;
+        SP2=Pt1;
+        NbPointsTotal=2;
       }
       else NbPointsTotal=3;///case SP1!=Pt1 && SP1!=Pt2!
     }
     else if( (NbPoints==2)&&(NbPointsTotal==2) ) {//there is also SP1!=SP2
       if( (SP1.CheckSameSP(Pt1))||(SP1.CheckSameSP(Pt2)) ) {
-       if( (SP2.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt2)) )
-         NbPointsTotal=2;
-       else NbPointsTotal=3;
+        if( (SP2.CheckSameSP(Pt1))||(SP2.CheckSameSP(Pt2)) )
+          NbPointsTotal=2;
+        else NbPointsTotal=3;
       }
       else NbPointsTotal=3;
     }
@@ -1571,113 +1354,11 @@ void TestNbPoints(const Standard_Integer ,
 }
 //=======================================================================
 //function : StartingPointsResearch
-//purpose  : 
-//=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch
-  (const Standard_Integer T1,
-   const Standard_Integer T2,
-   IntPolyh_StartPoint &SP1, 
-   IntPolyh_StartPoint &SP2) const 
-{
-  const IntPolyh_Point &P1=TPoints1[TTriangles1[T1].FirstPoint()];
-  const IntPolyh_Point &P2=TPoints1[TTriangles1[T1].SecondPoint()];
-  const IntPolyh_Point &P3=TPoints1[TTriangles1[T1].ThirdPoint()];
-  const IntPolyh_Point &Q1=TPoints2[TTriangles2[T2].FirstPoint()];
-  const IntPolyh_Point &Q2=TPoints2[TTriangles2[T2].SecondPoint()];
-  const IntPolyh_Point &Q3=TPoints2[TTriangles2[T2].ThirdPoint()];
-
-
-  /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
-     The sides are (e1,e2,e3) and (f1,f2,f3).
-     The normals are n1 and m1*/
-
-  const IntPolyh_Point  e1=P2-P1;
-  const IntPolyh_Point  e2=P3-P2;
-  const IntPolyh_Point  e3=P1-P3;
-
-  const IntPolyh_Point  f1=Q2-Q1;
-  const IntPolyh_Point  f2=Q3-Q2;
-  const IntPolyh_Point  f3=Q1-Q3;
-  
-
-  IntPolyh_Point nn1,mm1;
-  nn1.Cross(e1, e2); //normal of the first triangle
-  mm1.Cross(f1, f2); //normal of the  second triangle
-
-  Standard_Real nn1modulus, mm1modulus;
-  nn1modulus=sqrt(nn1.SquareModulus());
-  mm1modulus=sqrt(mm1.SquareModulus());
-
-  //-------------------------------------------------------
-  ///calculation of intersection points between two triangles
-  //-------------------------------------------------------
-  Standard_Integer NbPoints=0;
-  Standard_Integer NbPointsTotal=0;
-  IntPolyh_StartPoint Pt1,Pt2;
-
-
-    ///check T1 normal
-    if(Abs(nn1modulus)<MyConfusionPrecision){//10.0e-20) {
-
-    }
-    else {
-      const IntPolyh_Point n1=nn1.Divide(nn1modulus);
-      ///T2 edges with T1
-      if(NbPointsTotal<2) {
-       NbPoints=TriangleEdgeContact(1,1,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<2) {
-       NbPoints=TriangleEdgeContact(1,2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<2) {
-       NbPoints=TriangleEdgeContact(1,3,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-    }
-
-    ///check T2 normal
-    if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
-
-    }
-    else {
-      const IntPolyh_Point m1=mm1.Divide(mm1modulus);
-      ///T1 edges with T2  
-      if(NbPointsTotal<2) {
-       NbPoints=TriangleEdgeContact(2,1,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<2) { 
-       NbPoints=TriangleEdgeContact(2,2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<2) {
-       NbPoints=TriangleEdgeContact(2,3,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-    }
-
-  /*  if( (NbPointsTotal >1)&&( Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
-      &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) )*/
-  if( (NbPoints)&&(SP1.CheckSameSP(SP2)) )
-    NbPointsTotal=1;
-  
-  SP1.SetCoupleValue(T1,T2);
-  SP2.SetCoupleValue(T1,T2);
-  return (NbPointsTotal);
-}
-//=======================================================================
-//function : StartingPointsResearch2
 //purpose  : From  two  triangles compute intersection  points.
 //           If I found   more  than two intersection  points
 //           it  means that those triangle are coplanar
 //=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch2
+Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch
   (const Standard_Integer T1,
    const Standard_Integer T2,
    IntPolyh_StartPoint &SP1, 
@@ -1731,21 +1412,21 @@ Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch2
       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
       ///T2 edges with T1
       if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
     }
 
@@ -1757,21 +1438,21 @@ Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch2
       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
       ///T1 edges with T2
       if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if(NbPointsTotal<3) { 
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
     }
   //  no need already checked in  TestNbPoints
@@ -1793,139 +1474,11 @@ Standard_Integer IntPolyh_MaillageAffinage::StartingPointsResearch2
 }
 //=======================================================================
 //function : NextStartingPointsResearch
-//purpose  : 
-//=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch
-  (const Standard_Integer T1,
-   const Standard_Integer T2,
-   const IntPolyh_StartPoint &SPInit,
-   IntPolyh_StartPoint &SPNext) const 
-{
-  Standard_Integer NbPointsTotal=0;
-  if( (T1<0)||(T2<0) ) NbPointsTotal=0;
-  else {
-    const IntPolyh_Point &P1=TPoints1[TTriangles1[T1].FirstPoint()];
-    const IntPolyh_Point &P2=TPoints1[TTriangles1[T1].SecondPoint()];
-    const IntPolyh_Point &P3=TPoints1[TTriangles1[T1].ThirdPoint()];
-    const IntPolyh_Point &Q1=TPoints2[TTriangles2[T2].FirstPoint()];
-    const IntPolyh_Point &Q2=TPoints2[TTriangles2[T2].SecondPoint()];
-    const IntPolyh_Point &Q3=TPoints2[TTriangles2[T2].ThirdPoint()];
-    
-  /* The first triangle is (p1,p2,p3).  The other is (q1,q2,q3).
-     The sides are (e1,e2,e3) and (f1,f2,f3).
-     The normals are n1 and m1*/
-
-    const IntPolyh_Point  e1=P2-P1;
-    const IntPolyh_Point  e2=P3-P2;
-    const IntPolyh_Point  e3=P1-P3;
-    
-    const IntPolyh_Point  f1=Q2-Q1;
-    const IntPolyh_Point  f2=Q3-Q2;
-    const IntPolyh_Point  f3=Q1-Q3;
-    
-    IntPolyh_Point nn1,mm1;
-    nn1.Cross(e1, e2); //normal to the first triangle
-    mm1.Cross(f1, f2); //normal to the second triangle
-
-    Standard_Real nn1modulus, mm1modulus;
-    nn1modulus=sqrt(nn1.SquareModulus());
-    mm1modulus=sqrt(mm1.SquareModulus());
-     
-    //-------------------------------------------------
-    ///calculation of intersections points between triangles
-    //-------------------------------------------------
-    Standard_Integer NbPoints=0;
-    IntPolyh_StartPoint SP1,SP2;
-            
-    ///check T1 normal
-    if(Abs(nn1modulus)<MyConfusionPrecision) {//10.0e-20){
-
-    }
-    else {
-      const IntPolyh_Point n1=nn1.Divide(nn1modulus);
-      ///T2 edges with T1
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(1,1,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(1,2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(1,3,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-    }
-
-    ///check T2 normal
-    if(Abs(mm1modulus)<MyConfusionPrecision) {//10.0e-20){
-
-    }
-    else {
-      const IntPolyh_Point m1=mm1.Divide(mm1modulus);
-      ///T1 edges with T2
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(2,1,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(2,2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-      
-      if(NbPointsTotal<3) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact(2,3,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
-      }
-    }
-
-    if (NbPointsTotal==1) {
-      /*      if( (Abs(SP1.U1()-SPInit.U1())<MyConfusionPrecision)
-             &&(Abs(SP1.V1()-SPInit.V1())<MyConfusionPrecision) ) */
-      if(SP1.CheckSameSP(SP2))
-       NbPointsTotal=0;
-      else {
-
-       NbPointsTotal=0;
-      }
-    }
-
-    //    if ( (NbPointsTotal==2)&&( Abs(SP1.U1()-SPInit.U1())<MyConfusionPrecision)
-    //&&( Abs(SP1.V1()-SPInit.V1())<MyConfusionPrecision) ) {
-    if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SPInit)) ) {
-      NbPointsTotal=1;//SP1 et SPInit sont identiques
-      SPNext=SP2;
-    }
-    //    if( (NbPointsTotal==2)&&( Abs(SP2.U1()-SPInit.U1())<MyConfusionPrecision)
-    //&&( Abs(SP2.V1()-SPInit.V1())<MyConfusionPrecision) ) {
-    if( (NbPointsTotal==2)&&(SP2.CheckSameSP(SPInit)) ) {
-      NbPointsTotal=1;//SP2 et SPInit sont identiques
-      SPNext=SP1;
-    }
-    if(NbPointsTotal>1) {
-
-    }
-  }
-  SPNext.SetCoupleValue(T1,T2);
-  return (NbPointsTotal);
-}
-//=======================================================================
-//function : NextStartingPointsResearch2
 //purpose  : from  two triangles  and an intersection   point I
 //           seach the other point (if it exist).
 //           This function is used by StartPointChain
 //=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch2
+Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch
   (const Standard_Integer T1,
    const Standard_Integer T2,
    const IntPolyh_StartPoint &SPInit,
@@ -1982,21 +1535,21 @@ Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch2
       const IntPolyh_Point n1=nn1.Divide(nn1modulus);
       ///T2 edges with T1
       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.FirstEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,1,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q1,Q2,f1,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.SecondEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,2,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q2,Q3,f2,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if( (NbPointsTotal<3)&&(EdgeInit2!=Tri2.ThirdEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
-       TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(1,3,Tri1,Tri2,P1,P2,P3,e1,e2,e3,Q3,Q1,f3,n1,Pt1,Pt2);
+        TestNbPoints(1,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
     }
     ///check T2 normal
@@ -2007,29 +1560,29 @@ Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch2
       const IntPolyh_Point m1=mm1.Divide(mm1modulus);
       ///T1 edges with T2
       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.FirstEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,1,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P1,P2,e1,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.SecondEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,2,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P2,P3,e2,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
       
       if( (NbPointsTotal<3)&&(EdgeInit1!=Tri1.ThirdEdge()) ) {
-       IntPolyh_StartPoint Pt1,Pt2;
-       NbPoints=TriangleEdgeContact2(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
-       TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
+        IntPolyh_StartPoint Pt1,Pt2;
+        NbPoints=TriangleEdgeContact(2,3,Tri1,Tri2,Q1,Q2,Q3,f1,f2,f3,P3,P1,e3,m1,Pt1,Pt2);
+        TestNbPoints(2,NbPoints,NbPointsTotal,Pt1,Pt2,SP1,SP2);
       }
     }
       
     if (NbPointsTotal==1) {
       if(SP1.CheckSameSP(SPInit))
-       NbPointsTotal=0;
+        NbPointsTotal=0;
       else {
-       SPNext=SP1;
+        SPNext=SP1;
       }
     }
     else if( (NbPointsTotal==2)&&(SP1.CheckSameSP(SPInit)) ) {
@@ -2053,222 +1606,38 @@ Standard_Integer IntPolyh_MaillageAffinage::NextStartingPointsResearch2
 //purpose  : 
 //=======================================================================
 void CalculPtsInterTriEdgeCoplanaires(const Standard_Integer TriSurfID,
-                                     const IntPolyh_Point &NormaleTri,
-                                     const IntPolyh_Point &PE1,
-                                     const IntPolyh_Point &PE2,
-                                     const IntPolyh_Point &Edge,
-                                     const IntPolyh_Point &PT1,
-                                     const IntPolyh_Point &PT2,
-                                     const IntPolyh_Point &Cote,
-                                     const Standard_Integer CoteIndex,
-                                     IntPolyh_StartPoint &SP1,
-                                     IntPolyh_StartPoint &SP2,
-                                     Standard_Integer &NbPoints) 
+                                      const IntPolyh_Point &NormaleTri,
+                                      const IntPolyh_Triangle &Tri1,
+                                      const IntPolyh_Triangle &Tri2,
+                                      const IntPolyh_Point &PE1,
+                                      const IntPolyh_Point &PE2,
+                                      const IntPolyh_Point &Edge,
+                                      const Standard_Integer EdgeIndex,
+                                      const IntPolyh_Point &PT1,
+                                      const IntPolyh_Point &PT2,
+                                      const IntPolyh_Point &Cote,
+                                      const Standard_Integer CoteIndex,
+                                      IntPolyh_StartPoint &SP1,
+                                      IntPolyh_StartPoint &SP2,
+                                      Standard_Integer &NbPoints)
 {
-  IntPolyh_Point TestParalleles;
-  TestParalleles.Cross(Edge,Cote);
-  if(sqrt(TestParalleles.SquareModulus())<MyConfusionPrecision) {
-    IntPolyh_Point Per;
-    Per.Cross(NormaleTri,Cote);
-    Standard_Real p1p = Per.Dot(PE1);
-    Standard_Real p2p = Per.Dot(PE2);
-    Standard_Real p0p = Per.Dot(PT1);    
-    if ( ( (p1p>=p0p)&&(p2p<=p0p) )||( (p1p<=p0p)&&(p2p>=p0p) ) ) {
-      Standard_Real lambda=(p1p-p0p)/(p1p-p2p);
-      if (lambda<-MyConfusionPrecision) {
-
-      }
-      IntPolyh_Point PIE=PE1+Edge*lambda;      
-      Standard_Real alpha=RealLast();
-      if(Cote.X()!=0) alpha=(PIE.X()-PT1.X())/Cote.X();
-      else if (Cote.Y()!=0) alpha=(PIE.Y()-PT1.Y())/Cote.Y();
-      else if (Cote.Z()!=0) alpha=(PIE.Z()-PT1.Z())/Cote.Z();
-      else {
-
-      }
-      if (alpha<-MyConfusionPrecision) {
-
-      }
-      else {
-       if (NbPoints==0) {
-         SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
-         if (TriSurfID==1) {
-           SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-           SP1.SetUV2(PIE.U(),PIE.V());
-           SP1.SetEdge1(CoteIndex);
-           NbPoints++;
-         }
-         else if (TriSurfID==2) {
-           SP1.SetUV1(PIE.U(),PIE.V());
-           SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-           SP1.SetEdge2(CoteIndex);
-           NbPoints++;
-         }
-         else {
-
-         }
-       }
-
-       else if (NbPoints==1) {
-         SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
-         if (TriSurfID==1) {
-           SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-           SP2.SetUV2(PIE.U(),PIE.V());
-           SP2.SetEdge1(CoteIndex);
-           NbPoints++;
-         }
-         else if (TriSurfID==2) {
-           SP2.SetUV1(PIE.U(),PIE.V());
-           SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-           SP2.SetEdge2(CoteIndex);
-           NbPoints++;
-         }
-         else {
-
-         }
-       }
-
-       else if( (NbPoints>2)||(NbPoints<0) ) {
-
-       }
-      }
-    }
-  }
-  else {    //Cote et Edge paralleles, avec les rejections precedentes ils sont sur la meme droite
-    //On projette les points sur cette droite
-    Standard_Real pe1p= Cote.Dot(PE1);
-    Standard_Real pe2p= Cote.Dot(PE2);
-    Standard_Real pt1p= Cote.Dot(PT1);
-    Standard_Real pt2p= Cote.Dot(PT2);
-    
-    IntPolyh_Point PEP1,PTP1,PEP2,PTP2;
-
-    //PEP1 et PEP2 sont les points de contact entre le triangle et l'edge dans le repere UV de l'edge
-    //PTP1 et PTP2 sont les correspondants respectifs a PEP1 et PEP2 dans le repere UV du triangle
-
-
-    if (pe1p>pe2p) {
-      if ( (pt1p<pe1p)&&(pe1p<=pt2p) ) {
-       PEP1=PE1;
-       PTP1=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
-       NbPoints=1;
-       if (pt1p<=pe2p) {
-         PEP2=PE2;
-         PTP2=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
-         NbPoints=2;
-       }
-       else {
-         PEP2=PE1+Edge*((pt1p-pe1p)/(pe2p-pe1p));
-         PTP2=PT1;
-         NbPoints=2;
-       }
-      }
-      else if( (pt2p<pe1p)&&(pe1p<=pt1p) ) {
-       PEP1=PE1;
-       PTP1=PT1+Cote*((pt1p-pe1p)/(pt1p-pt2p));
-       NbPoints=1;
-       if (pt2p<=pe2p) {
-         PEP2=PE2;
-         PTP2=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
-         NbPoints=2;
-       }
-       else {
-         PEP2=PE1+Edge*((pt2p-pe1p)/(pe2p-pe1p));
-         PTP2=PT2;
-         NbPoints=2;
-       }
-      }
-    }
-    
-    if (pe1p<pe2p) {
-      if ( (pt1p<pe2p)&&(pe2p<=pt2p) ) {
-       PEP1=PE2;
-       PTP1=PT1+Cote*((pe2p-pt1p)/(pt2p-pt1p));
-       NbPoints=1;
-       if (pt1p<=pe1p) {
-         PEP2=PE1;
-         PTP2=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
-         NbPoints=2;
-       }
-       else {
-         PEP2=PE2+Edge*((pt1p-pe1p)/(pe2p-pe1p));
-         PTP2=PT1;
-         NbPoints=2;
-       }
-      }
-      else if( (pt2p<pe2p)&&(pe2p<=pt1p) ) {
-       PEP1=PE2;
-       PTP1=PT1+Cote*((pt1p-pe2p)/(pt1p-pt2p));
-       NbPoints=1;
-       if (pt2p<=pe1p) {
-         PEP2=PE1;
-         PTP2=PT1+Cote*((pe1p-pt1p)/(pt2p-pt1p));
-         NbPoints=2;
-       }
-       else {
-         PEP2=PE1+Edge*((pt2p-pe1p)/(pe2p-pe1p));
-         PTP2=PT2;
-         NbPoints=2;
-       }
-      }
-    }
-  
-    if (NbPoints!=0) {
-      if (Abs(PEP1.U()-PEP2.U())<MyConfusionPrecision
-          &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
-      
-      SP1.SetXYZ(PEP1.X(),PEP1.Y(),PEP1.Z());
-      if (TriSurfID==1) {
-       SP1.SetUV1(PTP1.U(),PTP1.V());    
-       SP1.SetUV2(PEP1.U(),PEP1.V());
-       SP1.SetEdge1(CoteIndex);
-      }
-      if (TriSurfID==2) {
-       SP1.SetUV1(PEP1.U(),PTP1.V());    
-       SP1.SetUV2(PTP1.U(),PEP1.V());
-       SP1.SetEdge2(CoteIndex);
-      }
-      
-      if (NbPoints==2) {
-       SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
-       if (TriSurfID==1) {
-         SP2.SetUV1(PTP2.U(),PTP2.V());          
-         SP2.SetUV2(PEP2.U(),PEP2.V());
-         SP2.SetEdge1(CoteIndex);
-       }
-       if (TriSurfID==2) {
-         SP2.SetUV1(PEP2.U(),PTP2.V());          
-         SP2.SetUV2(PTP2.U(),PEP2.V());
-         SP2.SetEdge2(CoteIndex);
-       } 
-      }
-    }
+  Standard_Real aDE, aDC;
+  //
+  gp_Vec aVE(Edge.X(), Edge.Y(), Edge.Z());
+  gp_Vec aVC(Cote.X(), Cote.Y(), Cote.Z());
+  //
+  aDE = aVE.SquareMagnitude();
+  aDC = aVC.SquareMagnitude();
+  //
+  if (aDE > SquareMyConfusionPrecision) {
+    aVE.Divide(aDE);
   }
-}
-//=======================================================================
-//function : CalculPtsInterTriEdgeCoplanaires2
-//purpose  : 
-//=======================================================================
-void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
-                                      const IntPolyh_Point &NormaleTri,
-                                      const IntPolyh_Triangle &Tri1,
-                                      const IntPolyh_Triangle &Tri2,
-                                      const IntPolyh_Point &PE1,
-                                      const IntPolyh_Point &PE2,
-                                      const IntPolyh_Point &Edge,
-                                      const Standard_Integer EdgeIndex,
-                                      const IntPolyh_Point &PT1,
-                                      const IntPolyh_Point &PT2,
-                                      const IntPolyh_Point &Cote,
-                                      const Standard_Integer CoteIndex,
-                                      IntPolyh_StartPoint &SP1,
-                                      IntPolyh_StartPoint &SP2,
-                                      Standard_Integer &NbPoints)
- {
-  IntPolyh_Point TestParalleles;
-  TestParalleles.Cross(Edge,Cote);
-
-  if(sqrt(TestParalleles.SquareModulus())>MyConfusionPrecision) {
+  //
+  if (aDC > SquareMyConfusionPrecision) {
+    aVC.Divide(aDC);
+  }
+  //
+  if (!aVE.IsParallel(aVC, MyConfusionPrecision)) {
     ///Edge and side are not parallel
     IntPolyh_Point Per;
     Per.Cross(NormaleTri,Cote);
@@ -2276,18 +1645,18 @@ void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
     Standard_Real p2p = Per.Dot(PE2);
     Standard_Real p0p = Per.Dot(PT1);
     ///The edge are PT1 are projected on the perpendicular of the side in the plane of the triangle
-    if ( ( (p1p>=p0p)&&(p0p>=p2p) )||( (p1p<=p0p)&&(p0p<=p2p) ) ) {
+    if ( ( ((p1p>=p0p)&&(p0p>=p2p) )||( (p1p<=p0p)&&(p0p<=p2p) )) && (Abs(p1p-p2p) > MyConfusionPrecision)) {
       Standard_Real lambda=(p1p-p0p)/(p1p-p2p);
       if (lambda<-MyConfusionPrecision) {
 
       }
       IntPolyh_Point PIE;
       if (Abs(lambda)<MyConfusionPrecision)//lambda=0
-       PIE=PE1;
+        PIE=PE1;
       else if (Abs(lambda)>1.0-MyConfusionPrecision)//lambda=1
-       PIE=PE2;
+        PIE=PE2;
       else
-       PIE=PE1+Edge*lambda;
+        PIE=PE1+Edge*lambda;
 
       Standard_Real alpha=RealLast();
       if(Cote.X()!=0) alpha=(PIE.X()-PT1.X())/Cote.X();
@@ -2301,103 +1670,103 @@ void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
 
       }
       else {
-       if (NbPoints==0) {
-         SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
-         if (TriSurfID==1) {
-           if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
-             SP1.SetUV1(PT1.U(),PT1.V());
-             SP1.SetUV1(PIE.U(),PIE.V());
-             SP1.SetEdge1(-1);
-           }
-           if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
-             SP1.SetUV1(PT2.U(),PT2.V());
-             SP1.SetUV1(PIE.U(),PIE.V());
-             SP1.SetEdge1(-1);
-           }
-           else {
-             SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-             SP1.SetUV2(PIE.U(),PIE.V());
-             SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
-             if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha);
-             else SP1.SetLambda1(1.0-alpha);
-           }
-           NbPoints++;
-         }
-         else if (TriSurfID==2) {
-           if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
-             SP1.SetUV1(PT1.U(),PT1.V());
-             SP1.SetUV1(PIE.U(),PIE.V());
-             SP1.SetEdge2(-1);
-           }
-           if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
-             SP1.SetUV1(PT2.U(),PT2.V());
-             SP1.SetUV1(PIE.U(),PIE.V());
-             SP1.SetEdge2(-1);
-           }
-           else {
-             SP1.SetUV1(PIE.U(),PIE.V());
-             SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-             SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
-             if (Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda2(alpha);
-             else SP1.SetLambda2(1.0-alpha);
-           }
-           NbPoints++;
-         }
-         else {
-
-         }
-       }
-
-       else if (NbPoints==1) {
-         SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
-         if (TriSurfID==1) {
-           if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
-             SP2.SetUV1(PT1.U(),PT1.V());
-             SP2.SetUV1(PIE.U(),PIE.V());
-             SP2.SetEdge1(-1);
-           }
-           if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
-             SP2.SetUV1(PT2.U(),PT2.V());
-             SP2.SetUV1(PIE.U(),PIE.V());
-             SP2.SetEdge1(-1);
-           }
-           else {
-             SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-             SP2.SetUV2(PIE.U(),PIE.V());
-             SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
-             if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha);
-             else SP2.SetLambda1(1.0-alpha);
-           }
-           NbPoints++;
-         }
-         else if (TriSurfID==2) {
-           if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
-             SP2.SetUV1(PT1.U(),PT1.V());
-             SP2.SetUV1(PIE.U(),PIE.V());
-             SP2.SetEdge2(-1);
-           }
-           if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
-             SP2.SetUV1(PT2.U(),PT2.V());
-             SP2.SetUV1(PIE.U(),PIE.V());
-             SP2.SetEdge2(-1);
-           }
-           else {
-             SP2.SetUV1(PIE.U(),PIE.V());
-             SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
-             SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
-             if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda2(alpha);
-             else SP2.SetLambda2(1.0-alpha);
-           }
-           NbPoints++;
-         }
-         else {
-
-         }
-       }
-
-       else if( (NbPoints>2)||(NbPoints<0) ) {
-
-         }
+        if (NbPoints==0) {
+          SP1.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
+          if (TriSurfID==1) {
+            if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
+              SP1.SetUV1(PT1.U(),PT1.V());
+              SP1.SetUV1(PIE.U(),PIE.V());
+              SP1.SetEdge1(-1);
+            }
+            if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
+              SP1.SetUV1(PT2.U(),PT2.V());
+              SP1.SetUV1(PIE.U(),PIE.V());
+              SP1.SetEdge1(-1);
+            }
+            else {
+              SP1.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
+              SP1.SetUV2(PIE.U(),PIE.V());
+              SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
+              if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha);
+              else SP1.SetLambda1(1.0-alpha);
+            }
+            NbPoints++;
+          }
+          else if (TriSurfID==2) {
+            if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
+              SP1.SetUV1(PT1.U(),PT1.V());
+              SP1.SetUV1(PIE.U(),PIE.V());
+              SP1.SetEdge2(-1);
+            }
+            if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
+              SP1.SetUV1(PT2.U(),PT2.V());
+              SP1.SetUV1(PIE.U(),PIE.V());
+              SP1.SetEdge2(-1);
+            }
+            else {
+              SP1.SetUV1(PIE.U(),PIE.V());
+              SP1.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
+              SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
+              if (Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda2(alpha);
+              else SP1.SetLambda2(1.0-alpha);
+            }
+            NbPoints++;
+          }
+          else {
+
+          }
+        }
+
+        else if (NbPoints==1) {
+          SP2.SetXYZ(PIE.X(),PIE.Y(),PIE.Z());
+          if (TriSurfID==1) {
+            if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
+              SP2.SetUV1(PT1.U(),PT1.V());
+              SP2.SetUV1(PIE.U(),PIE.V());
+              SP2.SetEdge1(-1);
+            }
+            if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
+              SP2.SetUV1(PT2.U(),PT2.V());
+              SP2.SetUV1(PIE.U(),PIE.V());
+              SP2.SetEdge1(-1);
+            }
+            else {
+              SP2.SetUV1(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
+              SP2.SetUV2(PIE.U(),PIE.V());
+              SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
+              if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha);
+              else SP2.SetLambda1(1.0-alpha);
+            }
+            NbPoints++;
+          }
+          else if (TriSurfID==2) {
+            if(Abs(alpha)<MyConfusionPrecision) {//alpha=0
+              SP2.SetUV1(PT1.U(),PT1.V());
+              SP2.SetUV1(PIE.U(),PIE.V());
+              SP2.SetEdge2(-1);
+            }
+            if(Abs(alpha)>1.0-MyConfusionPrecision) {//alpha=1
+              SP2.SetUV1(PT2.U(),PT2.V());
+              SP2.SetUV1(PIE.U(),PIE.V());
+              SP2.SetEdge2(-1);
+            }
+            else {
+              SP2.SetUV1(PIE.U(),PIE.V());
+              SP2.SetUV2(PT1.U()+Cote.U()*alpha,PT1.V()+Cote.V()*alpha);
+              SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
+              if (Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda2(alpha);
+              else SP2.SetLambda2(1.0-alpha);
+            }
+            NbPoints++;
+          }
+          else {
+
+          }
+        }
+
+        else if( (NbPoints>2)||(NbPoints<0) ) {
+
+          }
       }
     }
   }
@@ -2409,155 +1778,151 @@ void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
     Standard_Real pe2p= Cote.Dot(PE2);
     Standard_Real pt1p= Cote.Dot(PT1);
     Standard_Real pt2p= Cote.Dot(PT2);
-#ifndef DEB    
-    Standard_Real lambda1 =0.,lambda2 =0.,alpha1 =0.,alpha2 =0.;
-#else
-    Standard_Real lambda1,lambda2,alpha1,alpha2;
-#endif
+    Standard_Real lambda1=0., lambda2=0., alpha1=0., alpha2=0.;
     IntPolyh_Point PEP1,PTP1,PEP2,PTP2;
 
     if (pe1p>pe2p) {
       if ( (pt1p<pe1p)&&(pe1p<=pt2p) ) {
-       lambda1=0.0;
-       PEP1=PE1;
-       alpha1=((pe1p-pt1p)/(pt2p-pt1p));
-       PTP1=PT1+Cote*alpha1;
-       NbPoints=1;
-       if (pt1p<=pe2p) {
-         lambda2=1.0;
-         PEP2=PE2;
-         alpha2=((pe2p-pt1p)/(pt2p-pt1p));
-         PTP2=PT1+Cote*alpha2;
-         NbPoints=2;
-       }
-       else {
-         lambda2=((pt1p-pe1p)/(pe2p-pe1p));
-         PEP2=PE1+Edge*lambda2;
-         alpha2=0.0;
-         PTP2=PT1;
-         NbPoints=2;
-       }
+        lambda1=0.0;
+        PEP1=PE1;
+        alpha1=((pe1p-pt1p)/(pt2p-pt1p));
+        PTP1=PT1+Cote*alpha1;
+        NbPoints=1;
+        if (pt1p<=pe2p) {
+          lambda2=1.0;
+          PEP2=PE2;
+          alpha2=((pe2p-pt1p)/(pt2p-pt1p));
+          PTP2=PT1+Cote*alpha2;
+          NbPoints=2;
+        }
+        else {
+          lambda2=((pt1p-pe1p)/(pe2p-pe1p));
+          PEP2=PE1+Edge*lambda2;
+          alpha2=0.0;
+          PTP2=PT1;
+          NbPoints=2;
+        }
       }
       else if( (pt2p<pe1p)&&(pe1p<=pt1p) ) {
-       lambda1=0.0;
-       PEP1=PE1;
-       alpha1=((pt1p-pe1p)/(pt1p-pt2p));
-       PTP1=PT1+Cote*alpha1;
-       NbPoints=1;
-       if (pt2p<=pe2p) {
-         lambda2=1.0;
-         PEP2=PE2;
-         alpha2=((pe2p-pt1p)/(pt2p-pt1p));
-         PTP2=PT1+Cote*alpha2;
-         NbPoints=2;
-       }
-       else {
-         lambda2=((pt2p-pe1p)/(pe2p-pe1p));
-         PEP2=PE1+Edge*lambda2;
-         alpha2=1.0;
-         PTP2=PT2;
-         NbPoints=2;
-       }
+        lambda1=0.0;
+        PEP1=PE1;
+        alpha1=((pt1p-pe1p)/(pt1p-pt2p));
+        PTP1=PT1+Cote*alpha1;
+        NbPoints=1;
+        if (pt2p<=pe2p) {
+          lambda2=1.0;
+          PEP2=PE2;
+          alpha2=((pe2p-pt1p)/(pt2p-pt1p));
+          PTP2=PT1+Cote*alpha2;
+          NbPoints=2;
+        }
+        else {
+          lambda2=((pt2p-pe1p)/(pe2p-pe1p));
+          PEP2=PE1+Edge*lambda2;
+          alpha2=1.0;
+          PTP2=PT2;
+          NbPoints=2;
+        }
       }
     }
     
     if (pe1p<pe2p) {
       if ( (pt1p<pe2p)&&(pe2p<=pt2p) ) {
-       lambda1=1.0;
-       PEP1=PE2;
-       alpha1=((pe2p-pt1p)/(pt2p-pt1p));
-       PTP1=PT1+Cote*alpha1;
-       NbPoints=1;
-       if (pt1p<=pe1p) {
-         lambda2=0.0;
-         PEP2=PE1;
-         alpha2=((pe1p-pt1p)/(pt2p-pt1p));
-         PTP2=PT1+Cote*alpha2;
-         NbPoints=2;
-       }
-       else {
-         lambda2=((pt1p-pe1p)/(pe2p-pe1p));
-         PEP2=PE2+Edge*lambda2;
-         alpha2=0.0;
-         PTP2=PT1;
-         NbPoints=2;
-       }
+        lambda1=1.0;
+        PEP1=PE2;
+        alpha1=((pe2p-pt1p)/(pt2p-pt1p));
+        PTP1=PT1+Cote*alpha1;
+        NbPoints=1;
+        if (pt1p<=pe1p) {
+          lambda2=0.0;
+          PEP2=PE1;
+          alpha2=((pe1p-pt1p)/(pt2p-pt1p));
+          PTP2=PT1+Cote*alpha2;
+          NbPoints=2;
+        }
+        else {
+          lambda2=((pt1p-pe1p)/(pe2p-pe1p));
+          PEP2=PE2+Edge*lambda2;
+          alpha2=0.0;
+          PTP2=PT1;
+          NbPoints=2;
+        }
       }
       else if( (pt2p<pe2p)&&(pe2p<=pt1p) ) {
-       lambda1=1.0;
-       PEP1=PE2;
-       alpha1=((pt1p-pe2p)/(pt1p-pt2p));
-       PTP1=PT1+Cote*alpha1;
-       NbPoints=1;
-       if (pt2p<=pe1p) {
-         lambda2=0.0;
-         PEP2=PE1;
-         alpha2=((pe1p-pt1p)/(pt2p-pt1p));
-         PTP2=PT1+Cote*alpha2;
-         NbPoints=2;
-       }
-       else {
-         lambda2=((pt2p-pe1p)/(pe2p-pe1p));
-         PEP2=PE1+Edge*lambda2;
-         alpha2=1.0;
-         PTP2=PT2;
-         NbPoints=2;
-       }
+        lambda1=1.0;
+        PEP1=PE2;
+        alpha1=((pt1p-pe2p)/(pt1p-pt2p));
+        PTP1=PT1+Cote*alpha1;
+        NbPoints=1;
+        if (pt2p<=pe1p) {
+          lambda2=0.0;
+          PEP2=PE1;
+          alpha2=((pe1p-pt1p)/(pt2p-pt1p));
+          PTP2=PT1+Cote*alpha2;
+          NbPoints=2;
+        }
+        else {
+          lambda2=((pt2p-pe1p)/(pe2p-pe1p));
+          PEP2=PE1+Edge*lambda2;
+          alpha2=1.0;
+          PTP2=PT2;
+          NbPoints=2;
+        }
       }
     }
   
     if (NbPoints!=0) {
       SP1.SetXYZ(PEP1.X(),PEP1.Y(),PEP1.Z());
       if (TriSurfID==1) {///cote appartient a Tri1
-       SP1.SetUV1(PTP1.U(),PTP1.V());    
-       SP1.SetUV2(PEP1.U(),PEP1.V());
-       SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
+        SP1.SetUV1(PTP1.U(),PTP1.V());          
+        SP1.SetUV2(PEP1.U(),PEP1.V());
+        SP1.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
 
-       if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
-       else SP1.SetLambda1(1.0-alpha1);
-       
-       if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
-       else SP1.SetLambda2(1.0-lambda1);
+        if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
+        else SP1.SetLambda1(1.0-alpha1);
+        
+        if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
+        else SP1.SetLambda2(1.0-lambda1);
       }
       if (TriSurfID==2) {///cote appartient a Tri2
-       SP1.SetUV1(PEP1.U(),PTP1.V());    
-       SP1.SetUV2(PTP1.U(),PEP1.V());
-       SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
+        SP1.SetUV1(PEP1.U(),PTP1.V());          
+        SP1.SetUV2(PTP1.U(),PEP1.V());
+        SP1.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
 
-       if(Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
-       else SP1.SetLambda1(1.0-alpha1);
-       
-       if(Tri1.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
-       else SP1.SetLambda2(1.0-lambda1);
+        if(Tri2.GetEdgeOrientation(CoteIndex)>0) SP1.SetLambda1(alpha1);
+        else SP1.SetLambda1(1.0-alpha1);
+        
+        if(Tri1.GetEdgeOrientation(EdgeIndex)>0) SP1.SetLambda2(lambda1);
+        else SP1.SetLambda2(1.0-lambda1);
       }
       
       //It is checked if PEP1!=PEP2
       if ( (NbPoints==2)&&(Abs(PEP1.U()-PEP2.U())<MyConfusionPrecision)
-        &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
+         &&(Abs(PEP1.V()-PEP2.V())<MyConfusionPrecision) ) NbPoints=1;
       if (NbPoints==2) {
-       SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
-       if (TriSurfID==1) {
-         SP2.SetUV1(PTP2.U(),PTP2.V());          
-         SP2.SetUV2(PEP2.U(),PEP2.V());
-         SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
-
-         if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
-         else SP2.SetLambda1(1.0-alpha1);
-         
-         if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
-         else SP2.SetLambda2(1.0-lambda1);
-       }
-       if (TriSurfID==2) {
-         SP2.SetUV1(PEP2.U(),PTP2.V());          
-         SP2.SetUV2(PTP2.U(),PEP2.V());
-         SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
-
-         if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
-         else SP2.SetLambda1(1.0-alpha1);
-         
-         if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
-         else SP2.SetLambda2(1.0-lambda1);
-       } 
+        SP2.SetXYZ(PEP2.X(),PEP2.Y(),PEP2.Z());
+        if (TriSurfID==1) {
+          SP2.SetUV1(PTP2.U(),PTP2.V());          
+          SP2.SetUV2(PEP2.U(),PEP2.V());
+          SP2.SetEdge1(Tri1.GetEdgeNumber(CoteIndex));
+
+          if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
+          else SP2.SetLambda1(1.0-alpha1);
+          
+          if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
+          else SP2.SetLambda2(1.0-lambda1);
+        }
+        if (TriSurfID==2) {
+          SP2.SetUV1(PEP2.U(),PTP2.V());          
+          SP2.SetUV2(PTP2.U(),PEP2.V());
+          SP2.SetEdge2(Tri2.GetEdgeNumber(CoteIndex));
+
+          if(Tri1.GetEdgeOrientation(CoteIndex)>0) SP2.SetLambda1(alpha1);
+          else SP2.SetLambda1(1.0-alpha1);
+          
+          if(Tri2.GetEdgeOrientation(EdgeIndex)>0) SP2.SetLambda2(lambda1);
+          else SP2.SetLambda2(1.0-lambda1);
+        
       }
     }
   }
@@ -2583,299 +1948,12 @@ void CalculPtsInterTriEdgeCoplanaires2(const Standard_Integer TriSurfID,
       SP2.SetEdge2(-1);
   }
 }
+
 //=======================================================================
 //function : TriangleEdgeContact
 //purpose  : 
 //=======================================================================
 Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact
-  (const Standard_Integer TriSurfID,
-   const Standard_Integer EdgeIndex,
-   const IntPolyh_Point &PT1,
-   const IntPolyh_Point &PT2,
-   const IntPolyh_Point &PT3,
-   const IntPolyh_Point &Cote12,
-   const IntPolyh_Point &Cote23,
-   const IntPolyh_Point &Cote31,
-   const IntPolyh_Point &PE1,const IntPolyh_Point &PE2,
-   const IntPolyh_Point &Edge,
-   const IntPolyh_Point &NormaleT,
-   IntPolyh_StartPoint &SP1,IntPolyh_StartPoint &SP2) const 
-{
-
-  Standard_Real lambda =0.;
-  Standard_Real alpha =0.;
-  Standard_Real beta =0.;
-
-    
-  //The edge, on which the point is located, is known.
-  if (TriSurfID==1) {
-    SP1.SetEdge2(EdgeIndex);
-    SP2.SetEdge2(EdgeIndex);
-  }
-  else if (TriSurfID==2) {
-    SP1.SetEdge1(EdgeIndex);
-    SP2.SetEdge1(EdgeIndex);
-  }
-  else {
-
-  }
-
-/**The edge is projected on the normal of the triangle if yes 
-  --> free intersection (point I)--> start point is found*/
-  Standard_Integer NbPoints=0;
-
-  if(NormaleT.SquareModulus()==0) {
-
-  }
-  else if( (Cote12.SquareModulus()==0)
-       ||(Cote23.SquareModulus()==0)
-       ||(Cote31.SquareModulus()==0) ) {
-
-  }
-  else if(Edge.SquareModulus()==0) {
-
-  }
-  else {
-    Standard_Real pe1 = NormaleT.Dot(PE1);
-    Standard_Real pe2 = NormaleT.Dot(PE2);
-    Standard_Real pt1 = NormaleT.Dot(PT1);  
-    
-    // PE1I = lambda.Edge
-    
-    if( (Abs(pe1-pe2)<MyConfusionPrecision)&&(Abs(pe1-pt1)<MyConfusionPrecision) ) {
-      //edge and triangle are coplanar (two contact points maximum)
-
-      //The tops of the triangle are projected on the perpendicular of the edge 
-      
-      IntPolyh_Point PerpEdge;
-      PerpEdge.Cross(NormaleT,Edge);
-      Standard_Real pp1 = PerpEdge.Dot(PT1);
-      Standard_Real pp2 = PerpEdge.Dot(PT2);
-      Standard_Real pp3 = PerpEdge.Dot(PT3);
-      Standard_Real ppe1 = PerpEdge.Dot(PE1);
-               
-      if ( ( (pp1>ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
-       //there are two sides (commun top PT1) that can cut the edge
-       
-       //first side
-       CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                        PE1,PE2,Edge,PT1,PT2,
-                                        Cote12,1,SP1,SP2,NbPoints);
-
-       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-           &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-
-       //second side
-       if (NbPoints<2) {
-         CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                          PE1,PE2,Edge,PT3,PT1,
-                                          Cote31,3,SP1,SP2,NbPoints);
-       }
-      }
-
-      if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
-         &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-      if (NbPoints>=2) return(NbPoints);
-              
-      else if ( ( ( (pp2>ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
-              && (NbPoints<2) ) {
-       //there are two sides (common top PT2) that can cut the edge
-       
-       //first side
-       CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                        PE1,PE2,Edge,PT1,PT2,
-                                        Cote12,1,SP1,SP2,NbPoints);
-
-       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-           &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-
-       //second side
-       if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                                       PE1,PE2,Edge,PT2,PT3,
-                                                       Cote23,2,SP1,SP2,NbPoints);
-      }
-      if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-         &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-      if (NbPoints>=2) return(NbPoints);
-                     //= remove
-      else if ( ( (pp3>ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) )
-              && (NbPoints<2) ) {
-       //there are two sides (common top PT3) that can cut the edge
-       
-       //first side
-       CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                        PE1,PE2,Edge,PT3,PT1,Cote31,
-                                        3,SP1,SP2,NbPoints);
-       
-       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-           &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-
-       //second side
-       if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,
-                                                        PE1,PE2,Edge,PT2,PT3,
-                                                        Cote23,2,SP1,SP2,NbPoints);
-      }
-      if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-         &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
-      if (NbPoints>=2) return(NbPoints);
-    }
-    
-
-    //------------------------------------------------------
-    // edge and triangle NON COPLANAR (a contact point)
-    //------------------------------------------------------
-    else if(   ( (pe1>=pt1)&&(pe2<=pt1) ) || ( (pe1<=pt1)&&(pe2>=pt1) )  ) {
-      lambda=(pe1-pt1)/(pe1-pe2);
-      IntPolyh_Point PI;
-      if (lambda<-MyConfusionPrecision) {
-
-  }
-      else if (Abs(lambda)<MyConfusionPrecision) {//lambda==0
-       PI=PE1;
-       if(TriSurfID==1) SP1.SetEdge2(0);
-       else SP1.SetEdge1(0);
-      }
-      else if (Abs(lambda-1.0)<MyConfusionPrecision) {//lambda==1
-       PI=PE2;
-       if(TriSurfID==1) SP1.SetEdge2(0);
-       else SP1.SetEdge1(0);
-      }
-      else {
-       PI=PE1+Edge*lambda;
-       if(TriSurfID==1) SP1.SetEdge2(EdgeIndex);
-       else SP1.SetEdge1(EdgeIndex);
-      }
-     
-      if(Abs(Cote23.X())>MyConfusionPrecision) {
-       Standard_Real D=(Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23.X());
-       if(D!=0) alpha = (PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23.X())/D;
-       else { 
-
-       }
-       if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
-       else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
-      }
-      
-      else if (Abs(Cote12.X())>MyConfusionPrecision) { //On a Cote23.X()==0
-       alpha = (PI.X()-PT1.X())/Cote12.X();
-       
-       if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
-       
-       else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
-       else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
-       else  {
-
-       }
-      }
-      
-      else if (Abs(Cote23.Y())>MyConfusionPrecision) {
-       //On a Cote23.X()==0 et Cote12.X()==0 ==> equation can't be used
-       Standard_Real D=(Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y());
-       
-       if(D!=0) alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D;
-       else{ 
-
-       }
-  
-       if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
-       else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
-      }
-      
-      else if (Abs(Cote12.Y())>MyConfusionPrecision) {
-       //On a Cote23.X()==0, Cote12.X()==0 et Cote23.Y()==0
-       alpha = (PI.Y()-PT1.Y())/Cote12.Y();
-       
-       if ((Abs(alpha)<MyConfusionPrecision)||(Abs(alpha-1.0)<MyConfusionPrecision)) return(0);
-       
-       else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
-       
-       else {
-
-       }
-      }    
-      
-      else { //two equations of three can't be used
-
-       alpha=RealLast();
-       beta=RealLast();
-      }
-      
-      if( (beta<-MyConfusionPrecision)||(beta>(alpha+MyConfusionPrecision)) ) return(0);
-      else {
-       SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
-
-       if (TriSurfID==1) {
-         SP1.SetUV2(PI.U(),PI.V());
-         SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
-         NbPoints++;
-         if (beta<MyConfusionPrecision) {//beta==0 && alpha
-           SP1.SetEdge1(1);
-           SP1.SetLambda1(alpha);
-         }
-         if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
-           SP1.SetEdge1(3);
-           SP1.SetLambda1(1.0-alpha);
-         }
-         if (Abs(alpha-1)<MyConfusionPrecision)//alpha==1
-           SP1.SetEdge1(2);
-         if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
-           SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
-           SP1.SetUV1(PT1.U(),PT1.V());
-           SP1.SetEdge1(0);
-         }
-         if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
-           SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
-           SP1.SetUV1(PT2.U(),PT2.V());
-           SP1.SetEdge1(0);
-         }
-         if ( (Abs(beta-1)<MyConfusionPrecision)||(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
-           SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
-           SP1.SetUV1(PT3.U(),PT3.V());
-           SP1.SetEdge1(0);
-         }
-       }
-       else if(TriSurfID==2) {
-         SP1.SetUV1(PI.U(),PI.V());
-         SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
-         NbPoints++;
-         if (beta<MyConfusionPrecision) { //beta==0
-           SP1.SetEdge2(1);
-         }
-         if (Abs(beta-alpha)<MyConfusionPrecision)//beta==alpha
-           SP1.SetEdge2(3);
-         if (Abs(alpha-1)<MyConfusionPrecision)//alpha==1
-           SP1.SetEdge2(2);    
-         if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
-           SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
-           SP1.SetUV2(PT1.U(),PT1.V());
-           SP1.SetEdge2(0);
-         }
-         if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
-           SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
-           SP1.SetUV2(PT2.U(),PT2.V());
-           SP1.SetEdge2(0);
-         }
-         if ( (Abs(beta-1)<MyConfusionPrecision)||(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
-           SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
-           SP1.SetUV2(PT3.U(),PT3.V());
-           SP1.SetEdge2(0);
-         }
-       }
-       else{ 
-
-       }
-      }
-    }
-    else return 0;
-  }
-  return (NbPoints);
-}
-
-//=======================================================================
-//function : TriangleEdgeContact2
-//purpose  : 
-//=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
   (const Standard_Integer TriSurfID,
    const Standard_Integer EdgeIndex,
    const IntPolyh_Triangle &Tri1,
@@ -2945,62 +2023,62 @@ Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
 
       }
       else {
-       if ( ( (pp1>=ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<=ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
-         //there are two sides (common top PT1) that can cut the edge
-         
-         //first side
-         CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                           PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
-         
-         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-         
-         //second side
-         if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                                           PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
-       }
-       
-       if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
-           &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-       if (NbPoints>=2) return(NbPoints);
-       
-       else if ( ( ( (pp2>=ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<=ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
-                && (NbPoints<2) ) {
-         //there are two sides (common top PT2) that can cut the edge
-         
-         //first side
-         CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                           PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
-         
-         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-         
-         //second side
-         if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                                          PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
-       }
-       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-           &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-       if (NbPoints>=2) return(NbPoints);
-       
-       else if ( ( (pp3>=ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<=ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) )
-                && (NbPoints<2) ) {
-         //there are two sides (common top PT3) that can cut the edge
-         
-         //first side
-         CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                           PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
-         
-         if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-             &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
-         
-         //second side
-         if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires2(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
-                                                           PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
-       }
-       if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
-           &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
-       if (NbPoints>=2) return(NbPoints);
+        if ( ( (pp1>=ppe1)&&(pp2<=ppe1)&&(pp3<=ppe1) ) || ( (pp1<=ppe1)&&(pp2>=ppe1)&&(pp3>=ppe1) ) ){
+          //there are two sides (common top PT1) that can cut the edge
+          
+          //first side
+          CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                           PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
+          
+          if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
+              &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
+          
+          //second side
+          if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                                           PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
+        }
+        
+        if ( (NbPoints>1)&&(Abs(SP1.U1()-SP2.U1())<MyConfusionPrecision)
+            &&(Abs(SP1.V2()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
+        if (NbPoints>=2) return(NbPoints);
+        
+        else if ( ( ( (pp2>=ppe1)&&(pp1<=ppe1)&&(pp3<=ppe1) ) || ( (pp2<=ppe1)&&(pp1>=ppe1)&&(pp3>=ppe1) ) )
+                 && (NbPoints<2) ) {
+          //there are two sides (common top PT2) that can cut the edge
+          
+          //first side
+          CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                           PT1,PT2,Cote12,1,SP1,SP2,NbPoints);
+          
+          if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
+              &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
+          
+          //second side
+          if(NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                                          PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
+        }
+        if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
+            &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
+        if (NbPoints>=2) return(NbPoints);
+        
+        else if ( (( (pp3>=ppe1)&&(pp1<=ppe1)&&(pp2<=ppe1) ) || ( (pp3<=ppe1)&&(pp1>=ppe1)&&(pp2>=ppe1) ))
+                && (NbPoints<2) ) {
+          //there are two sides (common top PT3) that can cut the edge
+          
+          //first side
+          CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                           PT3,PT1,Cote31,3,SP1,SP2,NbPoints);
+          
+          if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
+              &&(Abs(SP1.V1()-SP2.V1())<MyConfusionPrecision) ) NbPoints=1;
+          
+          //second side
+          if (NbPoints<2) CalculPtsInterTriEdgeCoplanaires(TriSurfID,NormaleT,Tri1,Tri2,PE1,PE2,Edge,EdgeIndex,
+                                                           PT2,PT3,Cote23,2,SP1,SP2,NbPoints);
+        }
+        if ( (NbPoints>1)&&(Abs(SP2.U1()-SP1.U1())<MyConfusionPrecision)
+            &&(Abs(SP2.V1()-SP1.V1())<MyConfusionPrecision) ) NbPoints=1;
+        if (NbPoints>=2) return(NbPoints);
       }
     }
 
@@ -3014,26 +2092,27 @@ Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
 
       }
       else if (Abs(lambda)<MyConfusionPrecision) {//lambda==0
-       PI=PE1;
-       if(TriSurfID==1) SP1.SetEdge2(-1);
-       else SP1.SetEdge1(-1);
+        PI=PE1;
+        if(TriSurfID==1) SP1.SetEdge2(-1);
+        else SP1.SetEdge1(-1);
       }
       else if (Abs(lambda-1.0)<MyConfusionPrecision) {//lambda==1
-       PI=PE2;
-       if(TriSurfID==1) SP1.SetEdge2(-1);
-       else SP1.SetEdge1(-1);
+        PI=PE2;
+        if(TriSurfID==1) SP1.SetEdge2(-1);
+        else SP1.SetEdge1(-1);
       }
       else {
-       PI=PE1+Edge*lambda;
-       if(TriSurfID==1) 
-         if(Tri2.GetEdgeOrientation(EdgeIndex)>0)
-           SP1.SetLambda2(lambda);
-         else SP1.SetLambda2(1.0-lambda);
-       if(TriSurfID==2) 
-         if(Tri1.GetEdgeOrientation(EdgeIndex)>0)
-           SP1.SetLambda1(lambda);
-         else SP1.SetLambda1(1.0-lambda);
-
+        PI=PE1+Edge*lambda;
+        if(TriSurfID==1) {
+          if(Tri2.GetEdgeOrientation(EdgeIndex)>0)
+            SP1.SetLambda2(lambda);
+          else SP1.SetLambda2(1.0-lambda);
+        }
+        if(TriSurfID==2) {
+          if(Tri1.GetEdgeOrientation(EdgeIndex)>0)
+            SP1.SetLambda1(lambda);
+          else SP1.SetLambda1(1.0-lambda);
+        }
       }
         
       Standard_Real Cote23X=Cote23.X();
@@ -3042,164 +2121,164 @@ Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
 
       //Combination Eq1 Eq2
       if(Abs(Cote23X)>MyConfusionPrecision) {
-       D1=Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23X;
+        D1=Cote12.Y()-Cote12.X()*Cote23.Y()/Cote23X;
       }
       if(Abs(D1)>MyConfusionPrecision) {
-       alpha = ( PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23X )/D1;
+        alpha = ( PI.Y()-PT1.Y()-(PI.X()-PT1.X())*Cote23.Y()/Cote23X )/D1;
 
-       ///It is checked if 1.0>=alpha>=0.0
-       if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
-       else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23X;
+        ///It is checked if 1.0>=alpha>=0.0
+        if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
+        else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23X;
       }
       //Combination Eq1 and Eq2 with Cote23.X()==0
       else if ( (Abs(Cote12.X())>MyConfusionPrecision)
-             &&(Abs(Cote23X)<MyConfusionPrecision) ) { //There is Cote23.X()==0
-       alpha = (PI.X()-PT1.X())/Cote12.X();
-       
-       if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
-       
-       else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
-       else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
-       else  {
+              &&(Abs(Cote23X)<MyConfusionPrecision) ) { //There is Cote23.X()==0
+        alpha = (PI.X()-PT1.X())/Cote12.X();
+        
+        if ((alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision))) return(0);
+        
+        else if (Abs(Cote23.Y())>MyConfusionPrecision) beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
+        else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
+        else  {
 
-       }
+        }
       }
       //Combination Eq1 and Eq3
       else if ( (Abs(Cote23.X())>MyConfusionPrecision)
-             &&(Abs( D3= (Cote12.Z()-Cote12.X()*Cote23.Z()/Cote23.X()) ) > MyConfusionPrecision) ) {
-       
-       alpha = (PI.Z()-PT1.Z()-(PI.X()-PT1.X())*Cote23.Z()/Cote23.X())/D3;
-       
-       if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
-       else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
+              &&(Abs( D3= (Cote12.Z()-Cote12.X()*Cote23.Z()/Cote23.X()) ) > MyConfusionPrecision) ) {
+        
+        alpha = (PI.Z()-PT1.Z()-(PI.X()-PT1.X())*Cote23.Z()/Cote23.X())/D3;
+        
+        if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
+        else beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
       }
       //Combination Eq2 and Eq3
       else if ( (Abs(Cote23.Y())>MyConfusionPrecision)
-             &&(Abs( D4= (Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y()) ) > MyConfusionPrecision) ) {
-       
-       alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D4;
-       
-       if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
-       else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
+              &&(Abs( D4= (Cote12.Z()-Cote12.Y()*Cote23.Z()/Cote23.Y()) ) > MyConfusionPrecision) ) {
+        
+        alpha = (PI.Z()-PT1.Z()-(PI.Y()-PT1.Y())*Cote23.Z()/Cote23.Y())/D4;
+        
+        if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
+        else beta = (PI.Y()-PT1.Y()-alpha*Cote12.Y())/Cote23.Y();
       }
       //Combination Eq2 and Eq3 with Cote23.Y()==0
       else if ( (Abs(Cote12.Y())>MyConfusionPrecision)
-            && (Abs(Cote23.Y())<MyConfusionPrecision) ) {
-       alpha = (PI.Y()-PT1.Y())/Cote12.Y();
-       
-       if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
-       
-       else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
-       
-       else {
-         printf("\nCote PT2PT3 nul1\n");
-         PT2.Dump(2004);
-         PT3.Dump(3004);
-       }
+             && (Abs(Cote23.Y())<MyConfusionPrecision) ) {
+        alpha = (PI.Y()-PT1.Y())/Cote12.Y();
+        
+        if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
+        
+        else if (Abs(Cote23.Z())>MyConfusionPrecision) beta = (PI.Z()-PT1.Z()-alpha*Cote12.Z())/Cote23.Z();
+        
+        else {
+          printf("\nCote PT2PT3 nul1\n");
+          PT2.Dump(2004);
+          PT3.Dump(3004);
+        }
       }
       //Combination Eq1 and Eq3 with Cote23.Z()==0
       else if ( (Abs(Cote12.Z())>MyConfusionPrecision)
-            && (Abs(Cote23.Z())<MyConfusionPrecision) ) {
-       alpha = (PI.Z()-PT1.Z())/Cote12.Z();
-       
-       if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
-       
-       else if (Abs(Cote23.X())>MyConfusionPrecision) beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
-       
-       else {
+             && (Abs(Cote23.Z())<MyConfusionPrecision) ) {
+        alpha = (PI.Z()-PT1.Z())/Cote12.Z();
+        
+        if ( (alpha<-MyConfusionPrecision)||(alpha>(1.0+MyConfusionPrecision)) ) return(0);
+        
+        else if (Abs(Cote23.X())>MyConfusionPrecision) beta = (PI.X()-PT1.X()-alpha*Cote12.X())/Cote23.X();
+        
+        else {
 
-       }
+        }
       }
       
       else { //Particular case not processed ?
 
-       alpha=RealLast();
-       beta=RealLast();
+        alpha=RealLast();
+        beta=RealLast();
       }
       
       if( (beta<-MyConfusionPrecision)||(beta>(alpha+MyConfusionPrecision)) ) return(0);
       else {
-       SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
-
-       if (TriSurfID==1) {
-         SP1.SetUV2(PI.U(),PI.V());
-         SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
-         NbPoints++;
-         if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
-           SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
-           SP1.SetUV1(PT1.U(),PT1.V());
-           SP1.SetEdge1(-1);
-         }
-         else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
-           SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
-           SP1.SetUV1(PT2.U(),PT2.V());
-           SP1.SetEdge1(-1);
-         }
-         else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
-           SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
-           SP1.SetUV1(PT3.U(),PT3.V());
-           SP1.SetEdge1(-1);
-         }
-         else if (beta<MyConfusionPrecision) {//beta==0 
-           SP1.SetEdge1(Tri1.GetEdgeNumber(1));
-           if(Tri1.GetEdgeOrientation(1)>0)
-             SP1.SetLambda1(alpha);
-           else SP1.SetLambda1(1.0-alpha);
-         }
-         else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
-           SP1.SetEdge1(Tri1.GetEdgeNumber(3));
-           if(Tri1.GetEdgeOrientation(3)>0)
-             SP1.SetLambda1(1.0-alpha);
-           else SP1.SetLambda1(alpha);
-         }
-         else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
-           SP1.SetEdge1(Tri1.GetEdgeNumber(2));
-           if(Tri1.GetEdgeOrientation(2)>0)
-             SP1.SetLambda1(beta);
-           else SP1.SetLambda1(1.0-beta);
-         }
-       }
-       else if(TriSurfID==2) {
-         SP1.SetUV1(PI.U(),PI.V());
-         SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
-         NbPoints++;
-         if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
-           SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
-           SP1.SetUV2(PT1.U(),PT1.V());
-           SP1.SetEdge2(-1);
-         }
-         else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
-           SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
-           SP1.SetUV2(PT2.U(),PT2.V());
-           SP1.SetEdge2(-1);
-         }
-         else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
-           SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
-           SP1.SetUV2(PT3.U(),PT3.V());
-           SP1.SetEdge2(-1);
-         }
-         else if (beta<MyConfusionPrecision) { //beta==0
-           SP1.SetEdge2(Tri2.GetEdgeNumber(1));
-           if(Tri2.GetEdgeOrientation(1)>0)
-             SP1.SetLambda2(alpha);
-           else SP1.SetLambda2(1.0-alpha);
-         }
-         else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha
-           SP1.SetEdge2(Tri2.GetEdgeNumber(3));
-           if(Tri2.GetEdgeOrientation(3)>0)
-             SP1.SetLambda2(1.0-alpha);
-           else SP1.SetLambda2(alpha);
-         }
-         else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
-           SP1.SetEdge2(Tri2.GetEdgeNumber(2));        
-           if(Tri2.GetEdgeOrientation(2)>0)
-             SP1.SetLambda2(alpha);
-           else SP1.SetLambda2(1.0-alpha);
-         }
-       }
-       else {
-
-       }
+        SP1.SetXYZ(PI.X(),PI.Y(),PI.Z());
+
+        if (TriSurfID==1) {
+          SP1.SetUV2(PI.U(),PI.V());
+          SP1.SetUV1(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
+          NbPoints++;
+          if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
+            SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
+            SP1.SetUV1(PT1.U(),PT1.V());
+            SP1.SetEdge1(-1);
+          }
+          else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
+            SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
+            SP1.SetUV1(PT2.U(),PT2.V());
+            SP1.SetEdge1(-1);
+          }
+          else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
+            SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
+            SP1.SetUV1(PT3.U(),PT3.V());
+            SP1.SetEdge1(-1);
+          }
+          else if (beta<MyConfusionPrecision) {//beta==0 
+            SP1.SetEdge1(Tri1.GetEdgeNumber(1));
+            if(Tri1.GetEdgeOrientation(1)>0)
+              SP1.SetLambda1(alpha);
+            else SP1.SetLambda1(1.0-alpha);
+          }
+          else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha  
+            SP1.SetEdge1(Tri1.GetEdgeNumber(3));
+            if(Tri1.GetEdgeOrientation(3)>0)
+              SP1.SetLambda1(1.0-alpha);
+            else SP1.SetLambda1(alpha);
+          }
+          else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
+            SP1.SetEdge1(Tri1.GetEdgeNumber(2));
+            if(Tri1.GetEdgeOrientation(2)>0)
+              SP1.SetLambda1(beta);
+            else SP1.SetLambda1(1.0-beta);
+          }
+        }
+        else if(TriSurfID==2) {
+          SP1.SetUV1(PI.U(),PI.V());
+          SP1.SetUV2(PT1.U()+Cote12.U()*alpha+Cote23.U()*beta, PT1.V()+Cote12.V()*alpha+Cote23.V()*beta);
+          NbPoints++;
+          if (alpha<MyConfusionPrecision) {//alpha=0 --> beta==0
+            SP1.SetXYZ(PT1.X(),PT1.Y(),PT1.Z());
+            SP1.SetUV2(PT1.U(),PT1.V());
+            SP1.SetEdge2(-1);
+          }
+          else if ( (beta<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==0 alpha==1
+            SP1.SetXYZ(PT2.X(),PT2.Y(),PT2.Z());
+            SP1.SetUV2(PT2.U(),PT2.V());
+            SP1.SetEdge2(-1);
+          }
+          else if ( (Abs(beta-1)<MyConfusionPrecision)&&(Abs(1-alpha)<MyConfusionPrecision) ) {//beta==1 alpha==1
+            SP1.SetXYZ(PT3.X(),PT3.Y(),PT3.Z());
+            SP1.SetUV2(PT3.U(),PT3.V());
+            SP1.SetEdge2(-1);
+          }
+          else if (beta<MyConfusionPrecision) { //beta==0
+            SP1.SetEdge2(Tri2.GetEdgeNumber(1));
+            if(Tri2.GetEdgeOrientation(1)>0)
+              SP1.SetLambda2(alpha);
+            else SP1.SetLambda2(1.0-alpha);
+          }
+          else if (Abs(beta-alpha)<MyConfusionPrecision) {//beta==alpha
+            SP1.SetEdge2(Tri2.GetEdgeNumber(3));
+            if(Tri2.GetEdgeOrientation(3)>0)
+              SP1.SetLambda2(1.0-alpha);
+            else SP1.SetLambda2(alpha);
+          }
+          else if (Abs(alpha-1)<MyConfusionPrecision) {//alpha==1
+            SP1.SetEdge2(Tri2.GetEdgeNumber(2));        
+            if(Tri2.GetEdgeOrientation(2)>0)
+              SP1.SetLambda2(alpha);
+            else SP1.SetLambda2(1.0-alpha);
+          }
+        }
+        else {
+
+        }
       }
     }
     else return 0;
@@ -3207,72 +2286,6 @@ Standard_Integer IntPolyh_MaillageAffinage::TriangleEdgeContact2
   return (NbPoints);
 }
 //=======================================================================
-//function : TriangleComparePSP
-//purpose  : The   same as   TriangleCompare, plus compute the
-//           StartPoints without chaining them.
-//=======================================================================
-Standard_Integer IntPolyh_MaillageAffinage::TriangleComparePSP () 
-{
-  Standard_Integer CpteurTab=0;
-  Standard_Integer CpteurTabSP=0;
-  Standard_Real CoupleAngle=-2.0;
-  const Standard_Integer FinTT1 = TTriangles1.NbTriangles();
-  const Standard_Integer FinTT2 = TTriangles2.NbTriangles();
-
-  for(Standard_Integer i_S1=0; i_S1<FinTT1; i_S1++) {
-    for(Standard_Integer i_S2=0; i_S2<FinTT2; i_S2++){
-      if ( (TTriangles1[i_S1].IndiceIntersectionPossible() != 0)
-         &&(TTriangles1[i_S1].GetFleche() >= 0.0)
-         && (TTriangles2[i_S2].IndiceIntersectionPossible() != 0)
-         && (TTriangles2[i_S2].GetFleche() >= 0.0) ) {
-       IntPolyh_StartPoint SP1, SP2;
-       //If a triangle is dead or not in BSB, comparison is not possible
-       if (TriContact(TPoints1[TTriangles1[i_S1].FirstPoint()],
-                      TPoints1[TTriangles1[i_S1].SecondPoint()],
-                      TPoints1[TTriangles1[i_S1].ThirdPoint()],
-                      TPoints2[TTriangles2[i_S2].FirstPoint()],
-                      TPoints2[TTriangles2[i_S2].SecondPoint()],
-                      TPoints2[TTriangles2[i_S2].ThirdPoint()],
-                      CoupleAngle)){
-
-
-         TTriangles1[i_S1].SetIndiceIntersection(1);//The triangle is cut by another
-         TTriangles2[i_S2].SetIndiceIntersection(1);
-         
-         Standard_Integer NbPoints;
-         NbPoints=StartingPointsResearch(i_S1,i_S2,SP1, SP2);
-
-         if (NbPoints==0) {
-
-         }
-
-         if ( (NbPoints>0)&&(NbPoints<3) ) {
-           SP1.SetCoupleValue(i_S1,i_S2);
-           TStartPoints[CpteurTabSP]=SP1;
-           CpteurTabSP++;
-
-
-         }
-
-         if(NbPoints==2) {       
-           SP2.SetCoupleValue(i_S1,i_S2);
-           TStartPoints[CpteurTabSP]=SP2;
-           CpteurTabSP++;
-
-
-         }
-
-         if(NbPoints>2) {
-
-         }
-         CpteurTab++;
-       }
-      }
-    }
-  }
-  return(CpteurTabSP);
-}
-//=======================================================================
 //function : TriangleCompare
 //purpose  : Analyze  each couple of  triangles from the two --
 //           array  of triangles,  to   see  if they are  in
@@ -3281,151 +2294,121 @@ Standard_Integer IntPolyh_MaillageAffinage::TriangleComparePSP ()
 //=======================================================================
 Standard_Integer IntPolyh_MaillageAffinage::TriangleCompare ()
 {
-  Standard_Integer CpteurTab=0;
-
-  const Standard_Integer FinTT1 = TTriangles1.NbTriangles();
-  const Standard_Integer FinTT2 = TTriangles2.NbTriangles();
-
-  Standard_Integer TTClimit = 200;
-  Standard_Integer NbTTC = FinTT1 * FinTT2 / 10;
-  if (NbTTC < TTClimit)
-    NbTTC = TTClimit;
-  TTrianglesContacts.Init(NbTTC);
-  // eap
-  //TTrianglesContacts.Init(FinTT1 * FinTT2 / 10);
-
-  Standard_Real CoupleAngle=-2.0;
-  for(Standard_Integer i_S1=0; i_S1<FinTT1; i_S1++) {
-    for(Standard_Integer i_S2=0; i_S2<FinTT2; i_S2++){
-      if ( (TTriangles1[i_S1].IndiceIntersectionPossible() != 0)
-         &&(TTriangles1[i_S1].GetFleche() >= 0.0)
-         && (TTriangles2[i_S2].IndiceIntersectionPossible() != 0)
-         && (TTriangles2[i_S2].GetFleche() >= 0.0) ) {
-       //If a triangle is dead or not in BSB, comparison is not possible
-       IntPolyh_Triangle &Triangle1 =  TTriangles1[i_S1];
-       IntPolyh_Triangle &Triangle2 =  TTriangles2[i_S2];
-
-       if (TriContact(TPoints1[Triangle1.FirstPoint()],
-                      TPoints1[Triangle1.SecondPoint()],
-                      TPoints1[Triangle1.ThirdPoint()],
-                      TPoints2[Triangle2.FirstPoint()],
-                      TPoints2[Triangle2.SecondPoint()],
-                      TPoints2[Triangle2.ThirdPoint()],
-                      CoupleAngle)){
-
-         if (CpteurTab >= NbTTC)
-           {
-             TTrianglesContacts.SetNbCouples(CpteurTab);
-
-             return(CpteurTab);
-           }
-         TTrianglesContacts[CpteurTab].SetCoupleValue(i_S1, i_S2);
-         TTrianglesContacts[CpteurTab].SetAngleValue(CoupleAngle);
-//test   TTrianglesContacts[CpteurTab].Dump(CpteurTab);
-
-         Triangle1.SetIndiceIntersection(1);//The triangle is cut by another
-         Triangle2.SetIndiceIntersection(1);
-         CpteurTab++;
-       }
-      }
-    }
-  }
-  TTrianglesContacts.SetNbCouples(CpteurTab);
-
-  return(CpteurTab);
-}
-
-//=======================================================================
-//function : StartPointsCalcul
-//purpose  : From the array  of couples compute  all the start
-//           points and display them on the screen
-//=======================================================================
-void IntPolyh_MaillageAffinage::StartPointsCalcul() const
-{
-  const Standard_Integer FinTTC = TTrianglesContacts.NbCouples();
-//   printf("StartPointsCalcul() from IntPolyh_MaillageAffinage.cxx : StartPoints:\n");
-  for(Standard_Integer ii=0; ii<FinTTC; ii++) {
-    IntPolyh_StartPoint SP1,SP2;
-    Standard_Integer T1,T2;
-    T1=TTrianglesContacts[ii].FirstValue();
-    T2=TTrianglesContacts[ii].SecondValue();
-    StartingPointsResearch(T1,T2,SP1,SP2);
-    if ( (SP1.E1()!=-1)&&(SP1.E2()!=-1) ) SP1.Dump(ii);
-    if ( (SP2.E1()!=-1)&&(SP2.E2()!=-1) ) SP2.Dump(ii);
+  // Find couples with interfering bounding boxes
+  IntPolyh_IndexedDataMapOfIntegerArrayOfInteger aDMILI;
+  GetInterferingTriangles(TTriangles1, TPoints1,
+                          TTriangles2, TPoints2,
+                          aDMILI);
+  if (aDMILI.IsEmpty()) {
+    return 0;
+  }
+  //
+  Standard_Real CoupleAngle = -2.0;
+  //
+  // Intersection of the triangles
+  Standard_Integer i, aNb = aDMILI.Extent();
+  for (i = 1; i <= aNb; ++i) {
+    const Standard_Integer i_S1 = aDMILI.FindKey(i);
+    IntPolyh_Triangle &Triangle1 =  TTriangles1[i_S1];
+    const IntPolyh_Point& P1 = TPoints1[Triangle1.FirstPoint()];
+    const IntPolyh_Point& P2 = TPoints1[Triangle1.SecondPoint()];
+    const IntPolyh_Point& P3 = TPoints1[Triangle1.ThirdPoint()];
+    //
+    const IntPolyh_ArrayOfInteger& aLI2 = aDMILI(i);
+    IntPolyh_ArrayOfInteger::Iterator aItLI(aLI2);
+    for (; aItLI.More(); aItLI.Next()) {
+      const Standard_Integer i_S2 = aItLI.Value();
+      IntPolyh_Triangle &Triangle2 =  TTriangles2[i_S2];
+      const IntPolyh_Point& Q1 = TPoints2[Triangle2.FirstPoint()];
+      const IntPolyh_Point& Q2 = TPoints2[Triangle2.SecondPoint()];
+      const IntPolyh_Point& Q3 = TPoints2[Triangle2.ThirdPoint()];
+      //
+      if (TriContact(P1, P2, P3, Q1, Q2, Q3, CoupleAngle)) {
+        IntPolyh_Couple aCouple(i_S1, i_S2, CoupleAngle);
+        TTrianglesContacts.Append(aCouple);
+        //
+        Triangle1.SetIntersection(Standard_True);
+        Triangle2.SetIntersection(Standard_True);
+      }
+    }
   }
+  return TTrianglesContacts.Extent();
 }
+
 //=======================================================================
 //function : CheckCoupleAndGetAngle
 //purpose  : 
 //=======================================================================
-Standard_Integer CheckCoupleAndGetAngle(const Standard_Integer T1, 
-                                       const Standard_Integer T2,
-                                       Standard_Real& Angle, 
-                                       IntPolyh_ArrayOfCouples &TTrianglesContacts) 
+Standard_Boolean CheckCoupleAndGetAngle(const Standard_Integer T1, 
+                                        const Standard_Integer T2,
+                                        Standard_Real& Angle, 
+                                        IntPolyh_ListOfCouples &TTrianglesContacts) 
 {
-  Standard_Integer Test=0;
-  const Standard_Integer FinTTC = TTrianglesContacts.NbCouples();
-  for (Standard_Integer oioi=0; oioi<FinTTC; oioi++) {
-    IntPolyh_Couple TestCouple = TTrianglesContacts[oioi];
-    if ( (TestCouple.FirstValue()==T1)&&(TestCouple.AnalyseFlagValue()!=1) ) {
-      if (TestCouple.SecondValue()==T2) {
-       Test=oioi;
-       TTrianglesContacts[oioi].SetAnalyseFlag(1);
-       Angle=TTrianglesContacts[oioi].AngleValue();
-       oioi=FinTTC;
+  IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
+  for (; aIt.More(); aIt.Next()) {
+    IntPolyh_Couple& TestCouple = aIt.ChangeValue();
+    if (!TestCouple.IsAnalyzed()) {
+      if (TestCouple.FirstValue() == T1 && TestCouple.SecondValue() == T2) {
+        TestCouple.SetAnalyzed(Standard_True);
+        Angle = TestCouple.Angle();
+        return Standard_True;
       }
     }
   }
-  return(Test);
+  return Standard_False;
 }
 //=======================================================================
 //function : CheckCoupleAndGetAngle2
 //purpose  : 
 //=======================================================================
-Standard_Integer CheckCoupleAndGetAngle2(const Standard_Integer T1,
-                                        const Standard_Integer T2,
-                                        const Standard_Integer T11, 
-                                        const Standard_Integer T22,
-                                        Standard_Integer &CT11,
-                                        Standard_Integer &CT22, 
-                                        Standard_Real & Angle,
-                                        IntPolyh_ArrayOfCouples &TTrianglesContacts) 
+Standard_Boolean CheckCoupleAndGetAngle2(const Standard_Integer T1,
+                                         const Standard_Integer T2,
+                                         const Standard_Integer T11,
+                                         const Standard_Integer T22,
+                                         IntPolyh_ListIteratorOfListOfCouples& theItCT11,
+                                         IntPolyh_ListIteratorOfListOfCouples& theItCT22,
+                                         Standard_Real & Angle,
+                                         IntPolyh_ListOfCouples &TTrianglesContacts) 
 {
   ///couple T1 T2 is found in the list
   ///T11 and T22 are two other triangles implied  in the contact edge edge
   /// CT11 couple( T1,T22) and CT22 couple (T2,T11)
   /// these couples will be marked if there is a start point
-  Standard_Integer Test1=0;
-  Standard_Integer Test2=0;
-  Standard_Integer Test3=0;
-  const Standard_Integer FinTTC = TTrianglesContacts.NbCouples();
-  for (Standard_Integer oioi=0; oioi<FinTTC; oioi++) {
-    IntPolyh_Couple TestCouple = TTrianglesContacts[oioi];
-    if( (Test1==0)||(Test2==0)||(Test3==0) ) {
-      if ( (TestCouple.FirstValue()==T1)&&(TestCouple.AnalyseFlagValue()!=1) ) {
-       if (TestCouple.SecondValue()==T2) {
-         Test1=1;
-         TTrianglesContacts[oioi].SetAnalyseFlag(1);
-         Angle=TTrianglesContacts[oioi].AngleValue();
-       }
-       else if (TestCouple.SecondValue()==T22) {
-         Test2=1;
-         CT11=oioi;
-         Angle=TTrianglesContacts[oioi].AngleValue();
-       }
-      }
-      else if( (TestCouple.FirstValue()==T11)&&(TestCouple.AnalyseFlagValue()!=1) ) {
-       if (TestCouple.SecondValue()==T2) {
-         Test3=1;
-         CT22=oioi;
-         Angle=TTrianglesContacts[oioi].AngleValue();
-       }
-      }
-    }
-    else
-      oioi=FinTTC;
-  }
-  return(Test1);
+  Standard_Boolean Test1 , Test2, Test3;
+  Test1 = Test2 = Test3 = Standard_False;
+  //
+  IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
+  for (; aIt.More(); aIt.Next()) {
+    IntPolyh_Couple& TestCouple = aIt.ChangeValue();
+    if (TestCouple.IsAnalyzed()) {
+      continue;
+    }
+    //
+    if (TestCouple.FirstValue() == T1) {
+      if (TestCouple.SecondValue() == T2) {
+        Test1 = Standard_True;
+        TestCouple.SetAnalyzed(Standard_True);
+        Angle = TestCouple.Angle();
+      }
+      else if (TestCouple.SecondValue() == T22) {
+        Test2 = Standard_True;
+        theItCT11 = aIt;
+        Angle = TestCouple.Angle();
+      }
+    }
+    else if (TestCouple.FirstValue() == T11) {
+      if (TestCouple.SecondValue() == T2) {
+        Test3 = Standard_True;
+        theItCT22 = aIt;
+        Angle = TestCouple.Angle();
+      }
+    }
+    //
+    if (Test1 && Test2 && Test3) {
+      break;
+    }
+  }
+  return Test1;
 }
 //=======================================================================
 //function : CheckNextStartPoint
@@ -3434,30 +2417,30 @@ Standard_Integer CheckCoupleAndGetAngle2(const Standard_Integer T1,
 // the proper list number
 //=======================================================================
 Standard_Integer CheckNextStartPoint(IntPolyh_SectionLine & SectionLine,
-                                    IntPolyh_ArrayOfTangentZones & TTangentZones,
-                                    IntPolyh_StartPoint & SP,
-                                    const Standard_Boolean Prepend)//=Standard_False) 
+                                     IntPolyh_ArrayOfTangentZones & TTangentZones,
+                                     IntPolyh_StartPoint & SP,
+                                     const Standard_Boolean Prepend)//=Standard_False) 
 {
   Standard_Integer Test=1;
   if( (SP.E1()==-1)||(SP.E2()==-1) ) {
     //The tops of triangle are analyzed
     //It is checked if they are not in the array TTangentZones
-    Standard_Integer FinTTZ=TTangentZones.NbTangentZones();
+    Standard_Integer FinTTZ=TTangentZones.NbItems();
     for(Standard_Integer uiui=0; uiui<FinTTZ; uiui++) {
       IntPolyh_StartPoint TestSP=TTangentZones[uiui];
       if ( (Abs(SP.U1()-TestSP.U1())<MyConfusionPrecision)
-         &&(Abs(SP.V1()-TestSP.V1())<MyConfusionPrecision) ) {
-       if ( (Abs(SP.U2()-TestSP.U2())<MyConfusionPrecision)
-           &&(Abs(SP.V2()-TestSP.V2())<MyConfusionPrecision) ) {
-         Test=0;//SP is already in the list of  tops
-         uiui=FinTTZ;
-       }
+          &&(Abs(SP.V1()-TestSP.V1())<MyConfusionPrecision) ) {
+        if ( (Abs(SP.U2()-TestSP.U2())<MyConfusionPrecision)
+            &&(Abs(SP.V2()-TestSP.V2())<MyConfusionPrecision) ) {
+          Test=0;//SP is already in the list of  tops
+          uiui=FinTTZ;
+        }
       }
     }
     if (Test) {//the top does not belong to the list of TangentZones
       SP.SetChainList(-1);
       TTangentZones[FinTTZ]=SP;
-      TTangentZones.IncrementNbTangentZones();
+      TTangentZones.IncrementNbItems();
       Test=0;//the examined point is a top
     }
   }
@@ -3485,193 +2468,191 @@ Standard_Integer IntPolyh_MaillageAffinage::StartPointsChain
   (IntPolyh_ArrayOfSectionLines& TSectionLines,
    IntPolyh_ArrayOfTangentZones& TTangentZones) 
 {
-//Loop on the array of couples filled in the function COMPARE()
-  const Standard_Integer FinTTC = TTrianglesContacts.NbCouples();
-
-//Array of tops of triangles
-  for(Standard_Integer IndexA=0; IndexA<FinTTC; IndexA++) {
-    //First couple of triangles.
-    //It is checked if the couple of triangles has not been already examined.
-    if(TTrianglesContacts[IndexA].AnalyseFlagValue()!=1) {
-
-      Standard_Integer SectionLineIndex=TSectionLines.NbSectionLines();
+  //Loop on the array of couples filled in the function COMPARE()
+  IntPolyh_ListIteratorOfListOfCouples aIt(TTrianglesContacts);
+  for (; aIt.More(); aIt.Next()) {
+    IntPolyh_Couple& aCouple = aIt.ChangeValue();
+    // Check if the couple of triangles has not been already examined.
+    if(!aCouple.IsAnalyzed()) {
+
+      Standard_Integer SectionLineIndex=TSectionLines.NbItems();
       // fill last section line if still empty (eap)
       if (SectionLineIndex > 0
-         &&
-         TSectionLines[SectionLineIndex-1].NbStartPoints() == 0)
-       SectionLineIndex -= 1;
+          &&
+          TSectionLines[SectionLineIndex-1].NbStartPoints() == 0)
+        SectionLineIndex -= 1;
       else
-       TSectionLines.IncrementNbSectionLines();
+        TSectionLines.IncrementNbItems();
 
       IntPolyh_SectionLine &  MySectionLine=TSectionLines[SectionLineIndex];
       if (MySectionLine.GetN() == 0) // eap
-       MySectionLine.Init(10000);//Initialisation of array of StartPoint
+        MySectionLine.Init(10000);//Initialisation of array of StartPoint
 
       Standard_Integer NbPoints=-1;
       Standard_Integer T1I, T2I;
-      T1I = TTrianglesContacts[IndexA].FirstValue();
-      T2I = TTrianglesContacts[IndexA].SecondValue();
+      T1I = aCouple.FirstValue();
+      T2I = aCouple.SecondValue();
       
       // Start points for the current couple are found
       IntPolyh_StartPoint SP1, SP2;
-      NbPoints=StartingPointsResearch2(T1I,T2I,SP1, SP2);//first calculation
-      TTrianglesContacts[IndexA].SetAnalyseFlag(1);//the couple is marked
+      NbPoints=StartingPointsResearch(T1I,T2I,SP1, SP2);//first calculation
+      aCouple.SetAnalyzed(Standard_True);//the couple is marked
 
       if(NbPoints==1) {// particular case top/triangle or edge/edge
-       //the start point is input in the array
-       SP1.SetChainList(SectionLineIndex);
-       SP1.SetAngle(TTrianglesContacts[IndexA].AngleValue());
-       //it is checked if the point is not atop of the triangle
-       if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
-         IntPolyh_StartPoint SPNext1;
-         Standard_Integer TestSP1=0;
-         
-         //chain of a side
-         IntPolyh_StartPoint SP11;//=SP1;
-         if(SP1.E1()>=0) { //&&(SP1.E2()!=-1) already tested if the point is not a top
-           Standard_Integer NextTriangle1=0;
-           if (TEdges1[SP1.E1()].FirstTriangle()!=T1I) NextTriangle1=TEdges1[SP1.E1()].FirstTriangle();
-           else NextTriangle1=TEdges1[SP1.E1()].SecondTriangle();
-           
-           Standard_Real Angle=-2.0;
-           if (CheckCoupleAndGetAngle(NextTriangle1,T2I,Angle,TTrianglesContacts)) {
-             //it is checked if the couple exists and is marked
-             Standard_Integer NbPoints11=0;
-             NbPoints11=NextStartingPointsResearch2(NextTriangle1,T2I,SP1,SP11);
-             if (NbPoints11==1) {
-               SP11.SetChainList(SectionLineIndex);
-               SP11.SetAngle(Angle);
-               
-               if(CheckNextStartPoint(MySectionLine,TTangentZones,SP11)) {         
-                 Standard_Integer EndChainList=1;
-                 while (EndChainList!=0) {
-                   TestSP1=GetNextChainStartPoint(SP11,SPNext1,MySectionLine,TTangentZones);
-                   if(TestSP1==1) {
-                     SPNext1.SetChainList(SectionLineIndex);
-                     if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1))
-                       SP11=SPNext1;
-                     else EndChainList=0;
-                   }
-                   else EndChainList=0; //There is no next point
-                 }
-               }
-              
-             }
-             else {
-               if(NbPoints11>1) {//The point is input in the array TTangentZones
-                 TTangentZones[TTangentZones.NbTangentZones()]=SP11;//default list number = -1
-                 TTangentZones.IncrementNbTangentZones();
-               }
-               else {
-
-               }
-             }
-           }
-         }
-         else if (SP1.E2()<0){
-
-         }
-         //chain of the other side
-         IntPolyh_StartPoint SP12;//=SP1;
-         if (SP1.E2()>=0) { //&&(SP1.E1()!=-1) already tested
-           Standard_Integer NextTriangle2;
-           if (TEdges2[SP1.E2()].FirstTriangle()!=T2I) NextTriangle2=TEdges2[SP1.E2()].FirstTriangle();
-           else NextTriangle2=TEdges2[SP1.E2()].SecondTriangle();
-           
-           Standard_Real Angle=-2.0;
-           if(CheckCoupleAndGetAngle(T1I,NextTriangle2,Angle,TTrianglesContacts)) {
-             Standard_Integer NbPoints12=0;
-             NbPoints12=NextStartingPointsResearch2(T1I,NextTriangle2,SP1, SP12);
-             if (NbPoints12==1) {
-               
-               SP12.SetChainList(SectionLineIndex);
-               SP12.SetAngle(Angle);
-               Standard_Boolean Prepend = Standard_True; // eap
-               
-               if(CheckNextStartPoint(MySectionLine,TTangentZones,SP12, Prepend)) {
-                 Standard_Integer EndChainList=1;
-                 while (EndChainList!=0) {
-                   TestSP1=GetNextChainStartPoint(SP12,SPNext1,
-                                                  MySectionLine,TTangentZones,
-                                                  Prepend); // eap
-                   if(TestSP1==1) {
-                     SPNext1.SetChainList(SectionLineIndex);
-                     if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1,Prepend))
-                       SP12=SPNext1;
-                     else EndChainList=0;
-                   }
-                   else EndChainList=0; //there is no next point
-                 }
-               }
-               
-               else {
-                 if(NbPoints12>1) {//The points are input in the array TTangentZones
-                   TTangentZones[TTangentZones.NbTangentZones()]=SP12;//default list number = -1
-                   TTangentZones.IncrementNbTangentZones();
-                 }
-                 else {
-
-                 }
-               }
-             }
-           }
-         }
-         else if(SP1.E1()<0){
-
-         }
-       }
+        //the start point is input in the array
+        SP1.SetChainList(SectionLineIndex);
+        SP1.SetAngle(aCouple.Angle());
+        //it is checked if the point is not atop of the triangle
+        if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
+          IntPolyh_StartPoint SPNext1;
+          Standard_Integer TestSP1=0;
+          
+          //chain of a side
+          IntPolyh_StartPoint SP11;//=SP1;
+          if(SP1.E1()>=0) { //&&(SP1.E2()!=-1) already tested if the point is not a top
+            Standard_Integer NextTriangle1=0;
+            if (TEdges1[SP1.E1()].FirstTriangle()!=T1I) NextTriangle1=TEdges1[SP1.E1()].FirstTriangle();
+            else NextTriangle1=TEdges1[SP1.E1()].SecondTriangle();
+            
+            Standard_Real Angle=-2.0;
+            if (CheckCoupleAndGetAngle(NextTriangle1,T2I,Angle,TTrianglesContacts)) {
+              //it is checked if the couple exists and is marked
+              Standard_Integer NbPoints11=0;
+              NbPoints11=NextStartingPointsResearch(NextTriangle1,T2I,SP1,SP11);
+              if (NbPoints11==1) {
+                SP11.SetChainList(SectionLineIndex);
+                SP11.SetAngle(Angle);
+                
+                if(CheckNextStartPoint(MySectionLine,TTangentZones,SP11)) {            
+                  Standard_Integer EndChainList=1;
+                  while (EndChainList!=0) {
+                    TestSP1=GetNextChainStartPoint(SP11,SPNext1,MySectionLine,TTangentZones);
+                    if(TestSP1==1) {
+                      SPNext1.SetChainList(SectionLineIndex);
+                      if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1))
+                        SP11=SPNext1;
+                      else EndChainList=0;
+                    }
+                    else EndChainList=0; //There is no next point
+                  }
+                }
+               
+              }
+              else {
+                if(NbPoints11>1) {//The point is input in the array TTangentZones
+                  TTangentZones[TTangentZones.NbItems()]=SP11;//default list number = -1
+                  TTangentZones.IncrementNbItems();
+                }
+                else {
+
+                }
+              }
+            }
+          }
+          else if (SP1.E2()<0){
+
+          }
+          //chain of the other side
+          IntPolyh_StartPoint SP12;//=SP1;
+          if (SP1.E2()>=0) { //&&(SP1.E1()!=-1) already tested
+            Standard_Integer NextTriangle2;
+            if (TEdges2[SP1.E2()].FirstTriangle()!=T2I) NextTriangle2=TEdges2[SP1.E2()].FirstTriangle();
+            else NextTriangle2=TEdges2[SP1.E2()].SecondTriangle();
+            
+            Standard_Real Angle=-2.0;
+            if(CheckCoupleAndGetAngle(T1I,NextTriangle2,Angle,TTrianglesContacts)) {
+              Standard_Integer NbPoints12=0;
+              NbPoints12=NextStartingPointsResearch(T1I,NextTriangle2,SP1, SP12);
+              if (NbPoints12==1) {
+                
+                SP12.SetChainList(SectionLineIndex);
+                SP12.SetAngle(Angle);
+                Standard_Boolean Prepend = Standard_True; // eap
+                
+                if(CheckNextStartPoint(MySectionLine,TTangentZones,SP12, Prepend)) {
+                  Standard_Integer EndChainList=1;
+                  while (EndChainList!=0) {
+                    TestSP1=GetNextChainStartPoint(SP12,SPNext1,
+                                                   MySectionLine,TTangentZones,
+                                                   Prepend); // eap
+                    if(TestSP1==1) {
+                      SPNext1.SetChainList(SectionLineIndex);
+                      if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext1,Prepend))
+                        SP12=SPNext1;
+                      else EndChainList=0;
+                    }
+                    else EndChainList=0; //there is no next point
+                  }
+                }
+                
+                else {
+                  if(NbPoints12>1) {//The points are input in the array TTangentZones
+                    TTangentZones[TTangentZones.NbItems()]=SP12;//default list number = -1
+                    TTangentZones.IncrementNbItems();
+                  }
+                  else {
+
+                  }
+                }
+              }
+            }
+          }
+          else if(SP1.E1()<0){
+
+          }
+        }
       }
       else if(NbPoints==2) {
-       //the start points are input in the array 
-       IntPolyh_StartPoint SPNext2;
-       Standard_Integer TestSP2=0;
-       Standard_Integer EndChainList=1;
-
-       SP1.SetChainList(SectionLineIndex);
-       SP1.SetAngle(TTrianglesContacts[IndexA].AngleValue());
-       if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
-
-         //chain of a side
-         while (EndChainList!=0) {
-           TestSP2=GetNextChainStartPoint(SP1,SPNext2,MySectionLine,TTangentZones);
-           if(TestSP2==1) {
-             SPNext2.SetChainList(SectionLineIndex);
-             if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2))
-               SP1=SPNext2;
-             else EndChainList=0;
-           }
-           else EndChainList=0; //there is no next point
-         }
-       }
-       
+        //the start points are input in the array 
+        IntPolyh_StartPoint SPNext2;
+        Standard_Integer TestSP2=0;
+        Standard_Integer EndChainList=1;
+
+        SP1.SetChainList(SectionLineIndex);
+        SP1.SetAngle(aCouple.Angle());
+        if(CheckNextStartPoint(MySectionLine,TTangentZones,SP1)) {
+
+          //chain of a side
+          while (EndChainList!=0) {
+            TestSP2=GetNextChainStartPoint(SP1,SPNext2,MySectionLine,TTangentZones);
+            if(TestSP2==1) {
+              SPNext2.SetChainList(SectionLineIndex);
+              if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2))
+                SP1=SPNext2;
+              else EndChainList=0;
+            }
+            else EndChainList=0; //there is no next point
+          }
+        }
+        
         SP2.SetChainList(SectionLineIndex);
-       SP2.SetAngle(TTrianglesContacts[IndexA].AngleValue());
-       Standard_Boolean Prepend = Standard_True; // eap
-
-       if(CheckNextStartPoint(MySectionLine,TTangentZones,SP2,Prepend)) {
-         
-         //chain of the other side
-         EndChainList=1;
-         while (EndChainList!=0) {
-           TestSP2=GetNextChainStartPoint(SP2,SPNext2,
-                                          MySectionLine,TTangentZones,
-                                          Prepend); // eap
-           if(TestSP2==1) {
-             SPNext2.SetChainList(SectionLineIndex);
-             if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2,Prepend))
-               SP2=SPNext2;
-             else EndChainList=0;
-           }
-           else EndChainList=0; //there is no next point
-         }
-       }
+        SP2.SetAngle(aCouple.Angle());
+        Standard_Boolean Prepend = Standard_True; // eap
+
+        if(CheckNextStartPoint(MySectionLine,TTangentZones,SP2,Prepend)) {
+          
+          //chain of the other side
+          EndChainList=1;
+          while (EndChainList!=0) {
+            TestSP2=GetNextChainStartPoint(SP2,SPNext2,
+                                           MySectionLine,TTangentZones,
+                                           Prepend); // eap
+            if(TestSP2==1) {
+              SPNext2.SetChainList(SectionLineIndex);
+              if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext2,Prepend))
+                SP2=SPNext2;
+              else EndChainList=0;
+            }
+            else EndChainList=0; //there is no next point
+          }
+        }
       }
 
       else if( (NbPoints>2)&&(NbPoints<7) ) {
-       //More than two start points 
-       //the start point is input in the table
-       SP1.SetChainList(SectionLineIndex);
-       CheckNextStartPoint(MySectionLine,TTangentZones,SP1);
+        //More than two start points 
+        //the start point is input in the table
+        SP1.SetChainList(SectionLineIndex);
+        CheckNextStartPoint(MySectionLine,TTangentZones,SP1);
       }
       
       else {
@@ -3705,17 +2686,17 @@ Standard_Integer IntPolyh_MaillageAffinage::GetNextChainStartPoint
     //If is checked if two triangles intersect
     Standard_Real Angle= -2.0;
     if (CheckCoupleAndGetAngle(NextTriangle1,SP.T2(),Angle,TTrianglesContacts)) {
-      NbPoints=NextStartingPointsResearch2(NextTriangle1,SP.T2(),SP,SPNext);
+      NbPoints=NextStartingPointsResearch(NextTriangle1,SP.T2(),SP,SPNext);
       if( NbPoints!=1 ) {
-       if (NbPoints>1)
-         CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
-       else {
+        if (NbPoints>1)
+          CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
+        else {
 
-       NbPoints=0;
-       }
+        NbPoints=0;
+        }
       }
       else 
-       SPNext.SetAngle(Angle);
+        SPNext.SetAngle(Angle);
     }
     else NbPoints=0;//this couple does not intersect
   }
@@ -3727,17 +2708,17 @@ Standard_Integer IntPolyh_MaillageAffinage::GetNextChainStartPoint
       NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
     Standard_Real Angle= -2.0;
     if (CheckCoupleAndGetAngle(SP.T1(),NextTriangle2,Angle,TTrianglesContacts)) {
-      NbPoints=NextStartingPointsResearch2(SP.T1(),NextTriangle2,SP,SPNext);
+      NbPoints=NextStartingPointsResearch(SP.T1(),NextTriangle2,SP,SPNext);
       if( NbPoints!=1 ) {
-       if (NbPoints>1)
-         CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
-       else {
+        if (NbPoints>1)
+          CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend);
+        else {
 
-       NbPoints=0;
-       }
+        NbPoints=0;
+        }
       }
       else
-       SPNext.SetAngle(Angle);
+        SPNext.SetAngle(Angle);
     }
     else NbPoints=0;
   }
@@ -3749,43 +2730,36 @@ Standard_Integer IntPolyh_MaillageAffinage::GetNextChainStartPoint
   else if( (SP.E1()>=0)&&(SP.E2()>=0) ) {
     ///the point is located on two edges
       Standard_Integer NextTriangle1;
-      Standard_Integer CpleT11=-1;
-      Standard_Integer CpleT22=-1;
       if (TEdges1[SP.E1()].FirstTriangle()!=SP.T1()) NextTriangle1=TEdges1[SP.E1()].FirstTriangle();
       else 
-       NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
+        NextTriangle1=TEdges1[SP.E1()].SecondTriangle();
       Standard_Integer NextTriangle2;
       if (TEdges2[SP.E2()].FirstTriangle()!=SP.T2()) NextTriangle2=TEdges2[SP.E2()].FirstTriangle();
       else 
-       NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
+        NextTriangle2=TEdges2[SP.E2()].SecondTriangle();
       Standard_Real Angle= -2.0;
+
+      IntPolyh_ListIteratorOfListOfCouples aItCT11, aItCT22;
       if (CheckCoupleAndGetAngle2(NextTriangle1,NextTriangle2,
-                                 SP.T1(),SP.T2(),CpleT11,CpleT22,
-                                 Angle,TTrianglesContacts)) {
-       NbPoints=NextStartingPointsResearch2(NextTriangle1,NextTriangle2,SP,SPNext);
-       if( NbPoints!=1 ) {
-         if (NbPoints>1) {
-           ///The new point is checked
-           if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend)>0) {
-           }
-           else {
-
-           }
-         }
-         NbPoints=0;
-       }
-       else {//NbPoints==1
-         SPNext.SetAngle(Angle);     
-         //The couples (Ti,Tj) (Ti',Tj') are marked
-         if (CpleT11>=0) TTrianglesContacts[CpleT11].SetAnalyseFlag(1);
-         else {
-
-         }
-         if (CpleT22>=0) TTrianglesContacts[CpleT22].SetAnalyseFlag(1);
-         else {
-
-         }
-       }
+                                  SP.T1(),SP.T2(), aItCT11, aItCT22,
+                                  Angle,TTrianglesContacts)) {
+        NbPoints=NextStartingPointsResearch(NextTriangle1,NextTriangle2,SP,SPNext);
+        if( NbPoints!=1 ) {
+          if (NbPoints>1) {
+            ///The new point is checked
+            if(CheckNextStartPoint(MySectionLine,TTangentZones,SPNext,Prepend)>0) {
+            }
+            else {
+
+            }
+          }
+          NbPoints=0;
+        }
+        else {//NbPoints==1
+          SPNext.SetAngle(Angle);
+          if (aItCT11.More()) aItCT11.ChangeValue().SetAnalyzed(Standard_True);
+          if (aItCT22.More()) aItCT22.ChangeValue().SetAnalyzed(Standard_True);
+        }
       }
       else NbPoints=0;
     }
@@ -3839,26 +2813,11 @@ Bnd_Box IntPolyh_MaillageAffinage::GetBox(const Standard_Integer SurfID) const
   return(MyBox1);
   return(MyBox2);
 }
-
-//=======================================================================
-//function : GetBoxDraw
-//purpose  : 
-//=======================================================================
-void IntPolyh_MaillageAffinage::GetBoxDraw(const Standard_Integer SurfID)const
-{
-Standard_Real x0,y0,z0,x1,y1,z1;
-  if (SurfID==1) {
-    MyBox1.Get(x0,y0,z0,x1,y1,z1);
-  }
-  else {
-    MyBox2.Get(x0,y0,z0,x1,y1,z1);
-  }
-}
 //=======================================================================
 //function : GetArrayOfCouples
 //purpose  : 
 //=======================================================================
-IntPolyh_ArrayOfCouples &IntPolyh_MaillageAffinage::GetArrayOfCouples()
+IntPolyh_ListOfCouples &IntPolyh_MaillageAffinage::GetCouples()
 {
   return TTrianglesContacts;
 }
@@ -3866,7 +2825,7 @@ IntPolyh_ArrayOfCouples &IntPolyh_MaillageAffinage::GetArrayOfCouples()
 //function : SetEnlargeZone
 //purpose  : 
 //=======================================================================
-void IntPolyh_MaillageAffinage::SetEnlargeZone(Standard_Boolean& EnlargeZone)
+void IntPolyh_MaillageAffinage::SetEnlargeZone(const Standard_Boolean EnlargeZone)
 {
   myEnlargeZone = EnlargeZone;
 }
@@ -3878,17 +2837,35 @@ Standard_Boolean IntPolyh_MaillageAffinage::GetEnlargeZone() const
 {
   return myEnlargeZone;
 }
-//modified by NIZNHY-PKV Fri Jan 20 10:06:13 2012f
+
+//=======================================================================
+//function : GetMinDeflection
+//purpose  : 
+//=======================================================================
+Standard_Real IntPolyh_MaillageAffinage::GetMinDeflection(const Standard_Integer SurfID) const
+{
+  return (SurfID==1)? FlecheMin1:FlecheMin2;
+}
+
+//=======================================================================
+//function : GetMaxDeflection
+//purpose  : 
+//=======================================================================
+Standard_Real IntPolyh_MaillageAffinage::GetMaxDeflection(const Standard_Integer SurfID) const
+{
+  return (SurfID==1)? FlecheMax1:FlecheMax2;
+}
+
 //=======================================================================
 //function : DegeneratedIndex
 //purpose  : 
 //=======================================================================
 void DegeneratedIndex(const TColStd_Array1OfReal& aXpars,
-                     const Standard_Integer aNbX,
-                     const Handle(Adaptor3d_HSurface)& aS,
-                     const Standard_Integer aIsoDirection,
-                     Standard_Integer& aI1,
-                     Standard_Integer& aI2)
+                      const Standard_Integer aNbX,
+                      const Handle(Adaptor3d_HSurface)& aS,
+                      const Standard_Integer aIsoDirection,
+                      Standard_Integer& aI1,
+                      Standard_Integer& aI2)
 {
   Standard_Integer i;
   Standard_Boolean bDegX1, bDegX2;
@@ -3918,12 +2895,12 @@ void DegeneratedIndex(const TColStd_Array1OfReal& aXpars,
     aX=aXpars(i);
     if (bDegX1) {
       if (fabs(aX-aDegX1) < MyTolerance) {
-       aI1=i;
+        aI1=i;
       }
     }
     if (bDegX2) {
       if (fabs(aX-aDegX2) < MyTolerance) {
-       aI2=i;
+        aI2=i;
       }
     }
   }
@@ -3933,9 +2910,9 @@ void DegeneratedIndex(const TColStd_Array1OfReal& aXpars,
 //purpose  : 
 //=======================================================================
 Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
-                              const Standard_Integer aIndex,
-                              const Standard_Real aTol2,
-                              Standard_Real& aDegX)
+                               const Standard_Integer aIndex,
+                               const Standard_Real aTol2,
+                               Standard_Real& aDegX)
 {
   Standard_Boolean bRet;
   Standard_Integer i, aNbP;
@@ -3962,12 +2939,12 @@ Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
     for (i=1; i<aNbP; ++i) {
       aU=i*dU;
       if (i==aNbP-1){
-       aU=aU2;
+        aU=aU2;
       }
       aP2=aS->Value(aU, aV);
       aD2=aP1.SquareDistance(aP2);
       if (aD2>aTol2) {
-       return bRet;
+        return bRet;
       }
       aP1=aP2;
     }
@@ -3985,12 +2962,12 @@ Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
     for (i=1; i<aNbP; ++i) {
       aV=i*dV;
       if (i==aNbP-1){
-       aV=aV2;
+        aV=aV2;
       }
       aP2=aS->Value(aU, aV);
       aD2=aP1.SquareDistance(aP2);
       if (aD2>aTol2) {
-       return bRet;
+        return bRet;
       }
       aP1=aP2;
     }
@@ -4000,48 +2977,3 @@ Standard_Boolean IsDegenerated(const Handle(Adaptor3d_HSurface)& aS,
   //
   return bRet;
 }
-//modified by NIZNHY-PKV Fri Jan 20 10:06:15 2012t
-#ifdef DEB
-
-#include <TopoDS_Shape.hxx>
-#include <Poly_Triangulation.hxx>
-#include <TColgp_Array1OfPnt.hxx>
-#include <Poly_Array1OfTriangle.hxx>
-#include <BRep_TFace.hxx>
-#include <TopoDS_Face.hxx>
-
-//=======================================================================
-//function : TriangleShape
-//purpose  : shape with triangulation containing triangles
-//=======================================================================
-static TopoDS_Shape TriangleShape(const IntPolyh_ArrayOfTriangles & TTriangles,
-                                 const IntPolyh_ArrayOfPoints &    TPoints)
-{
-  TopoDS_Face aFace;
-  if (TPoints.NbPoints() < 1 || TTriangles.NbTriangles() < 1) return aFace;
-  
-  Handle(Poly_Triangulation) aPTriangulation =
-    new Poly_Triangulation(TPoints.NbPoints(),TTriangles.NbTriangles(),Standard_False);
-  TColgp_Array1OfPnt &       aPNodes         = aPTriangulation->ChangeNodes();
-  Poly_Array1OfTriangle &    aPTrialgles     = aPTriangulation->ChangeTriangles();
-  Standard_Integer i;
-  for (i=0; i<TPoints.NbPoints(); i++) {
-    const IntPolyh_Point& P = TPoints[i];
-    aPNodes(i+1).SetCoord(P.X(), P.Y(), P.Z());
-  }
-  for (i=0; i<TTriangles.NbTriangles(); i++) {
-    const IntPolyh_Triangle& T = TTriangles[i];
-    aPTrialgles(i+1).Set(T.FirstPoint()+1, T.SecondPoint()+1, T.ThirdPoint()+1);
-  }
-
-  Handle(BRep_TFace) aTFace = new BRep_TFace;
-  aTFace->Triangulation(aPTriangulation);
-  aFace.TShape(aTFace);
-  return aFace;
-}
-#endif
-
-//#define MyTolerance 10.0e-7
-//#define MyConfusionPrecision 10.0e-12
-//#define SquareMyConfusionPrecision 10.0e-24