--- ./OpenCASCADE6.3.0/ros/src/BRepMesh/BRepMesh_FastDiscret.cxx	2008-08-15 15:54:05.000000000 +0200
+++ ./oce/src/BRepMesh/BRepMesh_FastDiscret.cxx	2011-04-23 12:47:35.000000000 +0200
@@ -3,25 +3,21 @@
 // Author:	Ekaterina SMIRNOVA
 // Copyright: Open CASCADE SAS 2008
 
-
 #include <BRepMesh_FastDiscret.ixx>
 
-#include <BRepMesh_Array1OfVertexOfDelaun.hxx>
-#include <BRepMesh_ListIteratorOfListOfVertex.hxx>
+#include <BRepMesh_FastDiscretFace.hxx>
+#include <BRepMesh_FaceAttribute.hxx>
+#include <BRepMesh_DataStructureOfDelaun.hxx>
+#include <BRepMesh_ClassifierPtr.hxx>
 #include <BRepMesh_GeomTool.hxx>
-#include <MeshShape_ListIteratorOfListOfSurfacePoint.hxx>
-#include <MeshShape_CurvatureType.hxx>
-#include <MeshShape_Couple.hxx>
+#include <BRepMesh_PairOfPolygon.hxx>
+#include <BRepMesh_DataMapOfShapePairOfPolygon.hxx>
+#include <BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon.hxx>
 #include <GeomAbs_IsoType.hxx>
 #include <GeomAbs_SurfaceType.hxx>
 #include <TopAbs.hxx>
-//#include <TColStd_ListOfReal.hxx>
-#include <TColStd_Array1OfBoolean.hxx>
 #include <TColStd_HArray1OfReal.hxx>
-#include <TColStd_ListIteratorOfListOfReal.hxx>
-#include <TColStd_MapIteratorOfMapOfInteger.hxx>
 #include <Precision.hxx>
-#include <TColgp_Array1OfXY.hxx>
 
 #include <BRep_Builder.hxx>
 #include <BRep_Tool.hxx>
@@ -49,16 +45,12 @@
 
 #include <Geom2d_Curve.hxx>
 
-#include <MeshShape_DataMapOfShapeListOfTransient.hxx>
-#include <TColStd_ListOfTransient.hxx>
-#include <MeshShape_DataMapIteratorOfDataMapOfShapeListOfTransient.hxx>
 #include <TColStd_DataMapOfIntegerInteger.hxx>
 #include <BRepMesh_ShapeTool.hxx>
 #include <ElSLib.hxx>
 #include <Geom_Surface.hxx>
 #include <Adaptor3d_IsoCurve.hxx>
 #include <BRepMesh_IndexedMapOfVertex.hxx>
-#include <BRepMesh_Delaun.hxx>
 #include <Extrema_LocateExtPC.hxx>
 
 #include <BRepMesh_ListOfXY.hxx>
@@ -68,7 +60,7 @@
 #include <BRepMesh_IDMapOfNodeOfDataStructureOfDelaun.hxx>
 #include <Standard_ErrorHandler.hxx>
 #include <Standard_Failure.hxx>
-#include <MeshDS_DataMapOfInteger.hxx>
+//#include <TColStd_DataMapOfInteger.hxx>
 #include <TColGeom2d_SequenceOfCurve.hxx>
 #include <TopTools_SequenceOfShape.hxx>
 #include <NCollection_IncAllocator.hxx>
@@ -80,40 +72,20 @@
 #include <SortTools_ShellSortOfReal.hxx>
 #include <TCollection_CompareOfReal.hxx>
 
-static Standard_Integer LIMITE_TRIANGULATION = 1048576;
-
+#include <TopTools_HArray1OfShape.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
 
-static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
-{
-  Standard_Integer i, j, n = 0;
-  Standard_Real p, result = 0.;
+#include <vector>
 
-  for(i = 1; i <= PLen; i++)
-  {
-    // Sort
-    for(j = i + 1; j <= PLen; j++)
-    {
-      if(P(i) > P(j))
-      {
-        p = P(i);
-        P(i) = P(j);
-        P(j) = p;
-      }
-    }
-    // Accumulate
-    if (i != 1)
-    {
-      p = Abs(P(i) - P(i-1));
-      if(p > 1.e-7)
-      {
-        result += p;
-        n++;
-      }
-    }
-  }
-  return (n? (result / (Standard_Real) n) : -1.);
-}
-//--ofv.
+// NOTE: replaced by more correct check
+// #if defined(WNT) || defined(LIN)
+// #define HAVE_TBB 1
+// #endif
+
+// paralleling with Intel TBB
+#ifdef HAVE_TBB
+#include <tbb/parallel_for_each.h>
+#endif
 
 inline Standard_Real MaxFaceTol (const TopoDS_Face& theFace)
 {
@@ -141,13 +113,13 @@
 //purpose  : 
 //=======================================================================
 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    defle,
-                                   const Standard_Real    angl,
-                                   const Bnd_Box&         Box,
-				   const Standard_Boolean ws,
-				   const Standard_Boolean inshape,
-				   const Standard_Boolean relative,
-				   const Standard_Boolean shapetrigu): 
-  angle(angl), deflection(defle),WithShare(ws),nbDomains(0),  nbLocat(0), 
+                                           const Standard_Real    angl,
+                                           const Bnd_Box&         Box,
+                                           const Standard_Boolean ws,
+                                           const Standard_Boolean inshape,
+                                           const Standard_Boolean relative,
+                                           const Standard_Boolean shapetrigu): 
+  myangle(angl), mydeflection(defle),myWithShare(ws), nbLocat(0), 
   myrelative(relative), myshapetrigu(shapetrigu), myinshape(inshape),
   myInternalVerticesMode(Standard_True)
 {
@@ -164,21 +136,20 @@
   }
 }
 
-
 //=======================================================================
 //function : BRepMesh_FastDiscret
 //purpose  : 
 //=======================================================================
 
 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    defle,
-				   const TopoDS_Shape&    shape,
-                                   const Bnd_Box&         Box,
-				   const Standard_Real    angl,
-				   const Standard_Boolean ws,
-				   const Standard_Boolean inshape,
-				   const Standard_Boolean relative,
-				   const Standard_Boolean shapetrigu): 
- angle(angl), deflection(defle),WithShare(ws),nbDomains(0),  nbLocat(0), 
+                                           const TopoDS_Shape&    shape,
+                                           const Bnd_Box&         Box,
+                                           const Standard_Real    angl,
+                                           const Standard_Boolean ws,
+                                           const Standard_Boolean inshape,
+                                           const Standard_Boolean relative,
+                                           const Standard_Boolean shapetrigu): 
+ myangle(angl), mydeflection(defle),myWithShare(ws), nbLocat(0), 
  myrelative(relative), myshapetrigu(shapetrigu), myinshape(inshape),
  myInternalVerticesMode(Standard_True)
 {
@@ -193,25 +164,52 @@
     if (dy > mydtotale) mydtotale = dy;
     if (dz > mydtotale) mydtotale = dz;
   }
-  Add(shape);
+  Perform(shape);
 }
 
-
 //=======================================================================
-//function : Add(shape)
+//function : Perform(shape)
 //purpose  : 
 //=======================================================================
 
-void BRepMesh_FastDiscret::Add(const TopoDS_Shape& shape)
+void BRepMesh_FastDiscret::Perform(const TopoDS_Shape& shape)
 {
-  TopExp_Explorer ex;
-  for (ex.Init(shape, TopAbs_FACE); ex.More(); ex.Next()) {
-    Add(TopoDS::Face(ex.Current()));
+  std::vector<TopoDS_Face> aFaces;
+  for (TopExp_Explorer ex(shape, TopAbs_FACE); ex.More(); ex.Next()) {
+    TopoDS_Face aF = TopoDS::Face(ex.Current());
+    Add(aF);
+    aFaces.push_back(aF);
   }
+  
+  // mesh faces in parallel threads using TBB
+#ifdef HAVE_TBB
+  if (Standard::IsReentrant())
+    tbb::parallel_for_each (aFaces.begin(), aFaces.end(), *this);
+  else
+#endif
+  for (std::vector<TopoDS_Face>::iterator it(aFaces.begin()); it != aFaces.end(); it++)
+    Process (*it);
 }
 
 
 //=======================================================================
+//function : Process
+//purpose  : 
+//=======================================================================
+
+void BRepMesh_FastDiscret::Process(const TopoDS_Face& theFace) const
+{
+  //cout << "START face " << theFace.TShape().operator->() << endl << flush;
+  Handle(BRepMesh_FaceAttribute) fattribute;
+  if ( GetFaceAttribute (theFace, fattribute) ) 
+  {
+    BRepMesh_FastDiscretFace aTool (GetAngle(), WithShare(), InShape(), ShapeTrigu());
+    aTool.Add (theFace, fattribute, GetMapOfDefEdge());
+  }
+  //cout << "END   face " << theFace.TShape().operator->() << endl << flush;
+}
+
+//=======================================================================
 //function : Add(face)
 //purpose  : 
 //=======================================================================
@@ -233,13 +231,10 @@
   TopoDS_Face face = theface;
   BRepTools::Update(face);
   face.Orientation(TopAbs_FORWARD);
-  nbDomains++;
-  Domains.Bind(nbDomains, face);
   structure.Nullify();
   Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
   anAlloc->Reset(Standard_False);  
   structure=new BRepMesh_DataStructureOfDelaun(anAlloc);
-  structure->NewDomain(nbDomains);
   BRepAdaptor_Surface  BS(face, Standard_False);
   Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
   
@@ -250,36 +245,36 @@
 
   TopAbs_Orientation orFace = face.Orientation();
   Handle(Poly_Triangulation) T;
-  TopLoc_Location l;
+  TopLoc_Location loc;
 
-  if (!WithShare) {          
+  if (!myWithShare) {          
     vertices.Clear();
     edges.Clear();
   }
 
-  mylistver.Clear();
   myvemap.Clear();
   mylocation2d.Clear();
   internaledges.Clear();
 
-  Standard_Integer i, j;//,i1, i2, k;
+  Standard_Integer i;
   i = 1;
 
   Standard_Real defedge, defface;
   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, dx, dy, dz;
   Standard_Integer nbEdge = 0;
-  Standard_Real savangle = angle;
+  Standard_Real savangle = myangle;
   Standard_Real cdef;
   Standard_Real maxdef = 2.* MaxFaceTol(theface);
   defface = 0.;
 
-  if (!myrelative) defface = Max(deflection, maxdef);
+  if (!myrelative) defface = Max(mydeflection, maxdef);
 
   TColStd_SequenceOfReal aFSeq, aLSeq;
   TColGeom2d_SequenceOfCurve aCSeq;
   TopTools_SequenceOfShape aShSeq;
 
   TopoDS_Iterator exW(face);
+
   for (; exW.More(); exW.Next()) {
     const TopoDS_Shape& aWire = exW.Value();
     if (aWire.ShapeType() != TopAbs_WIRE)
@@ -291,9 +286,8 @@
       if (!mapdefle.IsBound(edge)) {
         if (myrelative) {
           if (edges.IsBound(edge)) {
-	    const TColStd_ListOfTransient& L = edges.Find(edge);
-	    const Handle(Poly_PolygonOnTriangulation)& P = 
-	      *(Handle(Poly_PolygonOnTriangulation)*)&(L.First());
+	    const BRepMesh_PairOfPolygon& pair = edges.Find(edge);
+	    const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
 	    defedge = P->Deflection();
 	  }
           else {
@@ -306,17 +300,17 @@
             defedge = dx;
             if (defedge < dy) defedge = dy;
             if (defedge < dz) defedge = dz;
-            // ajustement par rapport a la taille totale:
+            // adjusted in relation to the total size:
             cdef = mydtotale/(2*defedge);
             if (cdef < 0.5) cdef = 0.5;
             if (cdef > 2.) cdef = 2.;
-            defedge = cdef * defedge * deflection;
-            angle = savangle * cdef;
+            defedge = cdef * defedge * mydeflection;
+            myangle = savangle * cdef;
           }
           defface = defface + defedge;
  	  defface = Max(maxdef, defface);
         }
-        else defedge = deflection;
+        else defedge = mydeflection;
     
         defedge = Max(maxdef, defedge);
         defedge = Max(1.e-05 , defedge);
@@ -326,16 +320,16 @@
         defedge = mapdefle(edge);
         if (myrelative) {defface = defface + defedge; defface = Max(maxdef, defface);}
       }
-      Standard_Real f,l;
-      Handle(Geom2d_Curve) C = BRep_Tool::CurveOnSurface(edge, face, f, l);
+      Standard_Real f1,l1;
+      Handle(Geom2d_Curve) C = BRep_Tool::CurveOnSurface(edge, face, f1, l1);
       if (C.IsNull()) continue;
 
-      aFSeq.Append(f);
-      aLSeq.Append(l);
+      aFSeq.Append(f1);
+      aLSeq.Append(l1);
       aCSeq.Append(C);
       aShSeq.Append(edge);
-      Add(edge, face, gFace, C, defedge, f, l);
-      angle = savangle;
+      Add(edge, face, gFace, C, defedge, f1, l1);
+      myangle = savangle;
     }
   }
 
@@ -345,11 +339,11 @@
   }
 
   if (myrelative ) defface = defface / nbEdge;
-  else             defface = deflection;
+  else             defface = mydeflection;
 
-  if (WithShare) defface = Max(maxdef, defface);
+  if (myWithShare) defface = Max(maxdef, defface);
   
-  T = BRep_Tool::Triangulation(face, l);
+  T = BRep_Tool::Triangulation(face, loc);
 
   if (!myshapetrigu || T.IsNull()) {
     
@@ -359,9 +353,9 @@
     maxX=maxY=-1.e100;
     
     Standard_Integer ipn = 0;
-    Standard_Integer i =1;
-    for (i = 1; i <= myvemap.Extent(); i++) {
-      const BRepMesh_Vertex& aV = structure->GetNode(myvemap.FindKey(i));
+    Standard_Integer i1 =1;
+    for (i1 = 1; i1 <= myvemap.Extent(); i1++) {
+      const BRepMesh_Vertex& aV = structure->GetNode(myvemap.FindKey(i1));
       ipn++;
       xCur=aV.Coord().X();
       yCur=aV.Coord().Y();
@@ -370,10 +364,10 @@
       minY=Min(yCur, minY);
       maxY=Max(yCur, maxY);
     }
-    myumin = minX;
-    myumax = maxX;
-    myvmin = minY;
-    myvmax = maxY;
+    Standard_Real myumin = minX;
+    Standard_Real myumax = maxX;
+    Standard_Real myvmin = minY;
+    Standard_Real myvmax = maxY;
     
     const Standard_Real umin = BS.FirstUParameter();
     const Standard_Real umax = BS.LastUParameter();
@@ -412,9 +406,8 @@
       }
     }
 
-    // verif rapide de la validite des bornes calculees. Si fausse,
-    // surement un pb de pcurve.
-
+    // fast verification of the validity of calculated limits. If wrong, 
+    // sure a problem of pcurve.
     if (thetype == GeomAbs_BezierSurface &&
         (myumin < -0.5 || myumax > 1.5 || myvmin < -0.5 || myvmax > 1.5))
     {
@@ -423,14 +416,14 @@
  
     //define parameters for correct parametrics
     
-    deltaX = 1.0;
-    deltaY = 1.0;
+    Standard_Real deltaX = 1.0;
+    Standard_Real deltaY = 1.0;
     Standard_Integer nbVertices = myvemap.Extent();
     const Standard_Real tolclass = Precision::PConfusion(); //0.03*Max(myumax-myumin, myvmax-myvmin);
     
-    BRepMesh_Classifier* classifier = 
+    BRepMesh_ClassifierPtr classifier ( 
       new BRepMesh_Classifier(face, tolclass,	internaledges, myvemap, 
-                              structure, myumin, myumax, myvmin, myvmax);   
+                              structure, myumin, myumax, myvmin, myvmax) );   
   
     myfacestate = classifier->State();
     if (myfacestate == BRepMesh_SelfIntersectingWire)
@@ -442,23 +435,17 @@
         nbmaill++;
         
         //clear the structure of links
-        const MeshDS_MapOfInteger& aLinks = structure->LinkOfDomain(nbDomains);
-        for(MeshDS_MapOfInteger::Iterator anIter(aLinks);anIter.More();anIter.Next())
-          structure->ForseRemoveLink(anIter.Key());
-        for (i = 1; i <= myvemap.Extent(); i++)
-           structure->ForseRemoveNode(myvemap.FindKey(i));
+        structure.Nullify();
+        structure = new BRepMesh_DataStructureOfDelaun(anAlloc);
         
-        structure->ClearDeleted();
-        
-        mylistver.Clear();
         myvemap.Clear();
         mylocation2d.Clear();
         internaledges.Clear();
 
-        Standard_Integer j;
-        for(j = 1; j <= aShSeq.Length(); j++)
+        Standard_Integer j1;
+        for(j1 = 1; j1 <= aShSeq.Length(); j1++)
         {
-          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j));
+          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
           if (edges.IsBound(edge))
           {
             edges.UnBind(edge);
@@ -467,17 +454,17 @@
         }
         
         
-        for( j = 1; j <= aShSeq.Length(); j++)
+        for( j1 = 1; j1 <= aShSeq.Length(); j1++)
         {
-          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j));
+          const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
           defedge = mapdefle(edge) / 3.;
           defedge = Max(defedge, eps);
           mapdefle.Bind(edge, defedge);
-          const Handle(Geom2d_Curve)& C = aCSeq.Value(j);
-          Add(edge, face, gFace, C, defedge, aFSeq.Value(j), aLSeq.Value(j));
+          const Handle(Geom2d_Curve)& C = aCSeq.Value(j1);
+          Add(edge, face, gFace, C, defedge, aFSeq.Value(j1), aLSeq.Value(j1));
         }
 	  
-        delete classifier;
+        classifier.Nullify();
 
         classifier = new BRepMesh_Classifier(face, tolclass, internaledges, myvemap,
                                              structure, myumin, myumax, myvmin, myvmax);
@@ -493,24 +480,25 @@
     if (myfacestate != BRepMesh_NoError && myfacestate != BRepMesh_ReMesh)
     {
       mynottriangulated.Append(face);
-      delete classifier;
+      classifier.Nullify();
       return;
     }
     
+    //??
     //add internal vertices after self-intersection check
-    if(myInternalVerticesMode) {
+    /*if(myInternalVerticesMode) {
       for(TopExp_Explorer ex(face,TopAbs_VERTEX ,TopAbs_EDGE); ex.More(); ex.Next()) {
         const TopoDS_Vertex& aVert = TopoDS::Vertex(ex.Current());
         Add(aVert,face,gFace);
       }
       nbVertices = myvemap.Extent();
-    }
+    }*/
     
 
     
-    // essai de determination de la longueur vraie:
+    // try to find the real length:
     // akm (bug OCC16) : We must calculate these measures in non-singular
-    //     parts of face. Let`s try to compute average value of three
+    //     parts of face. Let's try to compute average value of three
     //     (umin, (umin+umax)/2, umax), and respectively for v.
     //                 vvvvv
     Standard_Real longu = 0.0, longv = 0.0; //, last , first;
@@ -524,7 +512,7 @@
     BS.D0 (myumin, myvmin, P11);
     BS.D0 (myumin, dfvave, P21);
     BS.D0 (myumin, myvmax, P31);
-    for (i=0, dfucur=myumin; i <= 20; i++, dfucur+=du)  {
+    for (i1=0, dfucur=myumin; i1 <= 20; i1++, dfucur+=du)  {
       BS.D0 (dfucur, myvmin, P12);
       BS.D0 (dfucur, dfvave, P22);
       BS.D0 (dfucur, myvmax, P32);
@@ -537,7 +525,7 @@
     BS.D0(myumin, myvmin, P11);
     BS.D0(dfuave, myvmin, P21);
     BS.D0(myumax, myvmin, P31);
-    for (i=0, dfvcur=myvmin; i <= 20; i++, dfvcur+=dv) {
+    for (i1=0, dfvcur=myvmin; i1 <= 20; i1++, dfvcur+=dv) {
       BS.D0 (myumin, dfvcur, P12);
       BS.D0 (dfuave, dfvcur, P22);
       BS.D0 (myumax, dfvcur, P32);
@@ -551,14 +539,12 @@
     // akm (bug OCC16) ^^^^^
     
     if (longu <= 1.e-16 || longv <= 1.e-16) {
-      //sisi, ca s`est vu!!
+      //yes, it is seen!!!
 #ifdef DEB_MESH_CHRONO
       chMaillEdges.Stop();
       MESH_CHRONO_TSTOP(thetype);
-      cout <<"WARNING: la face "<<myNbDomains<<" n'est pas triangulee."<<endl;
 #endif
-//    delete classifier;
-//    MESH_FAILURE(theface);
+      MESH_FAILURE(theface);
     }
 
 
@@ -569,7 +555,7 @@
       
       Dv = Max(1.0e0 - (defface/r),0.0e0) ;
       Standard_Real oldDv = 2.0 * ACos (Dv);
-      oldDv = Min(oldDv, angle);
+      oldDv = Min(oldDv, myangle);
       Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
       Dv = oldDv;
       
@@ -580,7 +566,7 @@
       if (ru > 1.e-16) {
         Du = Max(1.0e0 - (defface/ru),0.0e0);
         Du  = (2.0 * ACos (Du));
-        Du = Min(Du, angle);
+        Du = Min(Du, myangle);
         Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
         if(aa < gp::Resolution())
           return; 
@@ -593,294 +579,82 @@
       nbU = Max(nbU, (Standard_Integer)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
       
       Du = (myumax-myumin)/(nbU+1);
-      //-- On choisit DeltaX et DeltaY de facon a ce qu on ne saute pas 
-      //-- de points sur la grille
+      //-- DeltaX and DeltaY are chosen so that to avoid "jumping" 
+      //-- of points on the grid
       deltaX = Du;
       deltaY = Dv;
     }
     else if (thetype == GeomAbs_Cylinder) {
-      Standard_Real aMax = Max(myumax,myvmax);
+      /*Standard_Real aMax = Max(myumax,myvmax);
       deltaX = 0.01;  //1./aMax; //0.01;
-      deltaY = 1.0;
+      deltaY = 1.0;*/
+      gp_Cylinder Cyl = BS.Cylinder();
+      Standard_Real R = Cyl.Radius();
+      // Calculate parameters for iteration in U direction
+      Standard_Real Du = 1.0 - (defface/R);
+      if (Du < 0.0) Du = 0.0;
+      Du = 2.0 * ACos (Du);
+      if (Du > GetAngle()) Du = GetAngle();
+      deltaX = Du/longv;
+      deltaY = 1.;
     }
     else {
       deltaX = (myumax-myumin)/longu;
       deltaY = (myvmax-myvmin)/longv;
     }    
 
-
-    myUParam.Clear(); myVParam.Clear();
-        
-    TColStd_Array1OfInteger tabvert_corr(1, nbVertices);
-    
-    gp_Pnt2d p2d;
-    
-    // Check the necessity to fill the map of parameters
-    const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus ||
-                                         thetype == GeomAbs_BezierSurface ||
-                                         thetype == GeomAbs_BSplineSurface);
-    
-    myUParam.Clear(); myVParam.Clear();
-  
-    BRepMesh_IDMapOfNodeOfDataStructureOfDelaun aMoveNodes(myvemap.Extent());
-    //cout << "UVparam" << endl;
-    Standard_Real eps = Precision::PConfusion();
-    Standard_Integer nbu1 = 0, nbu2 = 0, nbv1 = 0, nbv2 = 0, nbtot = 0;
-    Standard_Boolean IsRectangle = Standard_False;
-    TColStd_ListOfReal U1Params, U2Params, V1Params, V2Params;
-
-    for (i = 1; i <= structure->NbNodes(); i++)
-    {
-      const BRepMesh_Vertex& v = structure->GetNode(i);
-      p2d = v.Coord();
-      if (useUVParam) {
-	//cout << p2d.X() << " " << p2d.Y() << endl;
-	Standard_Real u = p2d.X(), v = p2d.Y();
-        myUParam.Add(u);
-        myVParam.Add(v);
-
-	++nbtot;
-	if(Abs(u - myumin) <= eps) {
-	  ++nbu1;
-	  V1Params.Append(v);
-	}
-	if(Abs(u - myumax) <= eps) {
-	  ++nbu2;
-	  V2Params.Append(v);
-	}
-	if(Abs(v - myvmin) <= eps) {
-	  ++nbv1;
-	  U1Params.Append(u);
-	}
-	if(Abs(v - myvmax) <= eps) {
-	  ++nbv2;
-	  U2Params.Append(u);
-	}
-      }
-      gp_XY res;
-      res.SetCoord((p2d.X()-minX)/deltaX,(p2d.Y()-minY)/deltaY);
-      BRepMesh_Vertex v_new(res,v.Domain(),v.Location3d(),v.Movability());
-      const MeshDS_ListOfInteger& alist = structure->GetNodeList(i);
-      aMoveNodes.Add(v_new,alist);
-      tabvert_corr(i) = i;
-    }    
-    structure->ReplaceNodes(aMoveNodes);
-    
-    if(useUVParam) {
-      Standard_Integer nb = nbu1 + nbu2 + nbv1 + nbv2 - 4;
-      if((nbu1 == nbu2) && (nbv1 == nbv2) && (nbtot == nb)) IsRectangle = Standard_True;
-    }
-
-
-    Standard_Boolean rajout;
-
-     switch (thetype)
-     {
-     case GeomAbs_Plane:
-       rajout = !classifier->NaturalRestriction();
-       break;
-     case GeomAbs_Sphere:
-     case GeomAbs_Torus:
-       rajout = Standard_True;
-       break;
-     default:
-       rajout = Standard_False;
-     }  
-
-    BRepMesh_Delaun trigu( structure,tabvert_corr, orFace==TopAbs_FORWARD);
-
-    Standard_Boolean isaline;
-    isaline = ((myumax-myumin)<1.e-05) || ((myvmax-myvmin)<1.e-05);
-    
-    Standard_Real aDef = -1;
-    if (!isaline && structure->ElemOfDomain(nbDomains).Extent() > 0) {
-      TColStd_ListOfInteger badTri, nulTri;
-      
-      if(!rajout)
-      {
-        aDef = Control(gFace, defface, mylistver, badTri, nulTri, trigu, Standard_True);
-        if( aDef > defface || aDef < 0.)
-          rajout = Standard_True;
-      }
-
-      if(!rajout) {
-	if(useUVParam) {
-	  if(BS.IsUClosed()) {
-	    if(myVParam.Extent() > 2) {
-	      rajout = Standard_True;
-	    }
-	  }
-	  if(BS.IsVClosed()) {
-	    if(myUParam.Extent() > 2) {
-	      rajout = Standard_True;
-	    }
-	  }
-	}
-      }
-
-      if(rajout){
-	if(IsRectangle) {
-	  InternalVerticesForRectangle(gFace, mylistver, U1Params, U2Params, V1Params, V2Params);
-	}
-	else {
-	  InternalVertices(gFace, mylistver, defface,*classifier);
-	}
-
-        if (mylistver.Extent() > 0) {
-	  BRepMesh_Array1OfVertexOfDelaun verttab(1, mylistver.Extent());
-	  BRepMesh_ListIteratorOfListOfVertex itVer(mylistver);
-	  ipn = 1;
-	  for (; itVer.More(); itVer.Next())
-	    verttab(ipn++) = itVer.Value();
-	  trigu.AddVertices(verttab);
-	}
-        //control internal points
-        BRepMesh_ListOfVertex vvlist;
-        aDef = Control(gFace, defface, vvlist, badTri, nulTri, trigu, Standard_False);
-        mylistver.Append(vvlist);
-      }
+    if(!mapattrib.IsBound(theface))
+    {
+      Handle(BRepMesh_FaceAttribute) aFA = new BRepMesh_FaceAttribute();
+      aFA->GetDefFace() = defface;
+      aFA->GetUMax() = myumax;
+      aFA->GetVMax() = myvmax;
+      aFA->GetUMin() = myumin;
+      aFA->GetVMin() = myvmin;
+      aFA->GetDeltaX() = deltaX;
+      aFA->GetDeltaY() = deltaY;
+      aFA->GetMinX() = minX;
+      aFA->GetMinY() = minY;
+      aFA->GetClassifier() = classifier;
+      mapattrib.Bind(theface, aFA);
     }
 
-    //modify structure back
-    aMoveNodes.Clear();
-    for (i = 1; i <= structure->NbNodes(); i++)
-    {
-      const BRepMesh_Vertex& v = structure->GetNode(i);
-      p2d = v.Coord();
-      gp_XY res;
-      res.SetCoord(p2d.X()*deltaX+myumin,p2d.Y()*deltaY+myvmin);
-      BRepMesh_Vertex v_new(res,v.Domain(),v.Location3d(),v.Movability());
-      const MeshDS_ListOfInteger& alist = structure->GetNodeList(i);
-      aMoveNodes.Add(v_new,alist);
-    }
-    structure->ReplaceNodes(aMoveNodes);
+    //Standard_Integer nbVertices = myvemap.Extent();
+    T = new Poly_Triangulation(nbVertices, 1, Standard_True);
+    TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
+    TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
     
-    if (myinshape && T.IsNull()) AddInShape(face, (aDef < 0.0)? defface : aDef);
-  }
-
-  else {
-    // recup des triangles pour les mettre dans la structure.
-    const TColgp_Array1OfPnt& Nodes = T->Nodes();
-    const TColgp_Array1OfPnt2d& UVNodes = T->UVNodes();
-    const Poly_Array1OfTriangle& triangles = T->Triangles();
-    Standard_Integer nt, ed1, ed2, ed3, tri;
-#ifdef DEB
-    Standard_Integer nbNodes = 
-#endif
-      Nodes.Length();
-    Standard_Integer nbTriangles = T->NbTriangles();
-    BRepMesh_Vertex v1, v2, v3;
-    Standard_Integer iv1, iv2, iv3;
-    Standard_Integer k, niver, nFree;
-    Standard_Integer t[3], n[3];
-
-
-    Poly_Connect pc(T);
-
-    // count the free edges
-    nFree = 0;
-    TColStd_MapOfInteger vFree;
-    for (i = 1; i <= nbTriangles; i++) {
-      pc.Triangles(i,t[0],t[1],t[2]);
-      triangles(i).Get(n[0], n[1], n[2]);
-      for (j = 0; j < 3; j++) {
-	if (t[j] == 0)  {
-	  k = (j+1) % 3;
-	  if (!vFree.Contains(n[j])) vFree.Add(n[j]);
-	  if (!vFree.Contains(n[k])) vFree.Add(n[k]);
-	  nFree++;
-	}
-      }
-    }
-
-    TColStd_DataMapOfIntegerInteger IDMAP;
     Standard_Integer index;
-
-
-    // parmi ces noeuds exterieurs, trouver leur indice dans la structure.
-    TColStd_MapIteratorOfMapOfInteger it;
-
-    // iteration sur vemap;
-    for (i = 1; i <= myvemap.Extent(); i++) {
+    for (i = 1; i <= nbVertices; i++) {
       index = myvemap.FindKey(i);
+      Nodes(i) = Pnt(index);
+      Nodes2d(i).SetXY(Vertex(index).Coord());
+    }
 
-      for (it.Initialize(vFree); it.More(); it.Next()) {
-	const gp_Pnt2d& P2d = UVNodes(it.Key());
-	if (Vertex(index).Coord().IsEqual(P2d.Coord(), Precision::PConfusion())) {
-	  IDMAP.Bind(it.Key(), index);
-	  vFree.Remove(it.Key());
-	  break;
-	}
-      }
+    // storage of triangulation in the BRep.
+    //TopLoc_Location loc = face.Location();
+    if (!loc.IsIdentity()) {
+      gp_Trsf tr = loc.Transformation();
+      tr.Invert();
+      for (i = Nodes.Lower(); i <= Nodes.Upper(); i++) 
+	      Nodes(i).Transform(tr);
     }
-    
-    gp_Pnt P1, P2, P3;
 
-    for (nt = 1; nt <= nbTriangles; nt++) {
-      triangles(nt).Get(n[0], n[1], n[2]);
-      if (!l.IsIdentity()) {
-	P1 = Nodes(n[0]).Transformed(l.Transformation());
-	P2 = Nodes(n[1]).Transformed(l.Transformation());
-	P3 = Nodes(n[2]).Transformed(l.Transformation());
-      }
-      else {
-	P1 = Nodes(n[0]);
-	P2 = Nodes(n[1]);
-	P3 = Nodes(n[2]);
-      }
-      const gp_Pnt2d& P12d = UVNodes(n[0]);
-      const gp_Pnt2d& P22d = UVNodes(n[1]);
-      const gp_Pnt2d& P32d = UVNodes(n[2]);
-
-
-      if (!IDMAP.IsBound(n[0])) {
-	nbLocat++;
-	v1.Initialize(P12d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
-	Location3d.Bind(nbLocat, P1);
-	iv1 = structure->AddNode(v1);
-	IDMAP.Bind(n[0], iv1);
-      }
-      else {
-	iv1 = IDMAP.Find(n[0]);
-      }
-
-      if (!IDMAP.IsBound(n[1])) {
-	nbLocat++;
-	v2.Initialize(P22d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
-	Location3d.Bind(nbLocat, P2);
-	iv2 = structure->AddNode(v2);
-	IDMAP.Bind(n[1], iv2);
-      }
-      else {
-	iv2 = IDMAP.Find(n[1]);
-      }
-
-      if (!IDMAP.IsBound(n[2])) {
-	nbLocat++;
-	v3.Initialize(P32d.XY(), nbDomains, nbLocat, MeshDS_OnSurface);
-	Location3d.Bind(nbLocat, P3);
-	iv3 = structure->AddNode(v3);
-	IDMAP.Bind(n[2], iv3);
-      }
-      else {
-	iv3 = IDMAP.Find(n[2]);
-      }
-
-      if (face.Orientation() == TopAbs_REVERSED) {
-	niver=iv2;
-	iv2=iv3;
-	iv3=niver;
-      }
-
-      ed1=structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Free, nbDomains));
-      ed2=structure->AddLink(BRepMesh_Edge(iv2,iv3,MeshDS_Free, nbDomains));
-      ed3=structure->AddLink(BRepMesh_Edge(iv3,iv1,MeshDS_Free, nbDomains));
-      tri = structure->AddElement(BRepMesh_Triangle(Abs(ed1), Abs(ed2), Abs(ed3), 
-						    (ed1>0), (ed2>0), (ed3>0), 
-						    MeshDS_Free, nbDomains));
+    BRep_Builder B;
+    B.UpdateFace(face, T);
 
+    BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(internaledges);
+    for (; It.More(); It.Next()) {
+      const BRepMesh_PairOfPolygon& pair = It.Value();
+      const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
+      const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
+      if ( NOD1 == NOD2 )
+	B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
+      else
+	B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
     }
   }
+
 #ifndef DEB_MESH
   }
   catch(Standard_Failure)
@@ -888,6 +662,7 @@
     MESH_FAILURE(theface);
   }
 #endif // DEB_MESH
+  structure.Nullify();
 }
 
 
@@ -931,7 +706,7 @@
 
   const Standard_Boolean sameUV = uvFirst.IsEqual(uvLast, Precision::PConfusion());
 
-  //Controle vertice tolerances
+  //Control vertexes tolerances
   gp_Pnt pFirst = gFace->Value(uvFirst.X(), uvFirst.Y());
   gp_Pnt pLast = gFace->Value(uvLast.X(), uvLast.Y());
 
@@ -941,7 +716,7 @@
   if(mindist < BRep_Tool::Tolerance(pBegin) ||
      mindist < BRep_Tool::Tolerance(pEnd) ) mindist = defedge;
 
-  // controle des edges degenerees non codees.
+  // control of degenerated non-coded edges.
 
   Standard_Boolean degener = BRep_Tool::Degenerated(edge);
 
@@ -949,7 +724,7 @@
   {
     if (pBegin.IsSame(pEnd))
     {
-      // calcul de la longueur de l`edge en 3d
+      // calculation of the length of the edge in 3D
       Standard_Real longueur = 0.0;
       Standard_Real du = (wLast-wFirst)/20;
       gp_Pnt P1, P2;
@@ -974,7 +749,7 @@
   // Correct UV points
   if (sameUV)
   {
-    // 1. est-ce vraiment sameUV sans etre denegere
+    // 1. is it really sameUV without being degenerated
     gp_Pnt2d uvF, uvL;
     C2d->D0(first, uvF);
     C2d->D0(last, uvL);
@@ -1011,7 +786,7 @@
     vertices.Bind(pBegin, ipf);
   }
   theUV = FindUV(pBegin, uvFirst, ipf, gFace, mindist);
-  BRepMesh_Vertex vf(theUV, nbDomains, ipf, MeshDS_Frontier);
+  BRepMesh_Vertex vf(theUV, ipf, MeshDS_Frontier);
   Standard_Integer ivf = structure->AddNode(vf);
 
   // Process last vertex
@@ -1042,7 +817,7 @@
     }
   }
   theUV = FindUV(pEnd, uvLast, ipl, gFace, mindist);
-  BRepMesh_Vertex vl(theUV, nbDomains, ipl, MeshDS_Frontier);
+  BRepMesh_Vertex vl(theUV, ipl, MeshDS_Frontier);
   Standard_Integer ivl= structure->AddNode(vl);
 
   Standard_Integer isvf = myvemap.FindIndex(ivf);
@@ -1063,7 +838,7 @@
 
     if (degener)
     {
-      // creation du nouveau:
+      // creation anew:
       TColStd_Array1OfInteger Nodes(1, 2), NodInStruct(1, 2);
       TColStd_Array1OfReal Param(1, 2);
       
@@ -1085,60 +860,101 @@
       BRepAdaptor_Curve cons;
       if (BRep_Tool::SameParameter(edge))
       {
-	cons.Initialize(edge);
+	      cons.Initialize(edge);
       }
       else
       {
-	cons.Initialize(edge, face);
+	      cons.Initialize(edge, face);
       }
 
       TopLoc_Location L;
       Standard_Integer nbpmin = 2;
       if (cons.GetType() == GeomAbs_Circle) nbpmin = 4; //OCC287
-      BRepMesh_GeomTool GT(cons, wFirst, wLast, 0.5*angle, otherdefedge, nbpmin);
-      // Creation des polygones sur triangulation:
-      const Standard_Integer nbnodes = GT.NbPoints();
+      BRepMesh_GeomTool GT(cons, wFirst, wLast, 0.5*myangle, otherdefedge, nbpmin);
+
+      // PTv, chl/922/G9, Take into account internal vertices
+      // it is necessary for internal edges, which does not split other edges, by their vertex
+      TopoDS_Iterator exV(edge);
+      for ( ; exV.More(); exV.Next() )
+      {
+	      TopoDS_Vertex aIntV = TopoDS::Vertex(exV.Value());
+	      if ( aIntV.Orientation() == TopAbs_INTERNAL )
+         GT.AddPoint(BRep_Tool::Pnt(aIntV),
+	  	      BRep_Tool::Parameter(aIntV, edge),
+		        Standard_True);
+      }
+
+      // Creation of polygons on triangulation:
+      Standard_Real puv;
+      Standard_Integer i;      
+      Standard_Real aEdgeLen = 0., a2dEdgeLen = 0.;
+      gp_Pnt P3dPrev;
+      gp_Pnt2d uvPrev;
+      GT.Value(cons, gFace, 1, puv, P3dPrev, uvPrev);
+      for (i = 2; i <= GT.NbPoints(); i++)
+      {
+        // Record 3d point
+	      GT.Value(cons, gFace, i, puv, P3d, uv);
+        aEdgeLen += P3d.SquareDistance(P3dPrev);
+        a2dEdgeLen += uv.SquareDistance(uvPrev);
+        P3dPrev = P3d;
+        uvPrev = uv;
+      }
+
+      Standard_Integer nbnodes = GT.NbPoints();
+      Standard_Real aMinToler = Min(defedge, mydeflection);
+      if(aEdgeLen < aMinToler*aMinToler)
+      {
+        //check 2D resolution
+        Standard_Real ddu = gFace->UResolution(aMinToler);
+        Standard_Real ddv = gFace->VResolution(aMinToler);
+        Standard_Real aMin = Min(ddu, ddv);
+
+        if(a2dEdgeLen < aMin*aMin)
+          nbnodes = 2;
+      }
       TColStd_Array1OfInteger Nodes(1, nbnodes);
       TColStd_Array1OfInteger NodInStruct(1, nbnodes);
       TColStd_Array1OfReal Param(1, nbnodes);
-      
-      // traitement du 1er point:
+
+      // processing of the 1st point:
       Nodes(1) = isvf;
       NodInStruct(1) = ipf;
       Param(1) = wFirst;
       
-      // dernier point:
+      // last point:
       Nodes(nbnodes) = isvl;
       NodInStruct(nbnodes) = ipl;
       Param(nbnodes) = wLast;
 
-      Standard_Real puv;
-      Standard_Integer i;
-      Standard_Integer iv2;
-      for (i = 2; i < GT.NbPoints(); i++)
+      if(nbnodes > 2)
       {
-        // Record 3d point
-	GT.Value(cons, gFace, i, puv, P3d, uv);
-	nbLocat++;
-	Location3d.Bind(nbLocat, P3d);
-	NodInStruct(i) = nbLocat;
-        // Record 2d point
-	v2.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_OnCurve);
-        iv2=structure->AddNode(v2);
-        
-	Standard_Integer isv = myvemap.FindIndex(iv2);
-	if (isv == 0) isv = myvemap.Add(iv2);
-	Nodes(i) = isv;
-	NodInStruct(i) = nbLocat;
-	Param(i) = puv;
-	
-	if (orEdge == TopAbs_FORWARD)
-	  structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
-	else if (orEdge == TopAbs_REVERSED)
-	  structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
-	else if (orEdge == TopAbs_INTERNAL)
-	  structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
-	ivf = iv2;
+        Standard_Integer iv2;
+        for (i = 2; i < GT.NbPoints(); i++)
+        {
+          // Record 3d point
+	        GT.Value(cons, gFace, i, puv, P3d, uv);
+        	nbLocat++;
+        	Location3d.Bind(nbLocat, P3d);
+        	NodInStruct(i) = nbLocat;
+          // Record 2d point
+	        v2.Initialize(uv.Coord(), nbLocat, MeshDS_OnCurve);
+          iv2=structure->AddNode(v2);
+          
+          Standard_Integer isv = myvemap.FindIndex(iv2);
+          if (isv == 0) isv = myvemap.Add(iv2);
+          Nodes(i) = isv;
+          NodInStruct(i) = nbLocat;
+          Param(i) = puv;
+
+          if (orEdge == TopAbs_FORWARD)
+            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier));
+          else if (orEdge == TopAbs_REVERSED)
+            structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier));
+          else if (orEdge == TopAbs_INTERNAL)
+            structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed));
+          ivf = iv2;
+        }
       }
 
       P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
@@ -1146,23 +962,23 @@
     }
 
     P2->Deflection(otherdefedge);
-    TColStd_ListOfTransient L;
-    L.Append(P2);
-    edges.Bind(edge, L);
+    BRepMesh_PairOfPolygon pair;
+    pair.Append(P2);
+    edges.Bind(edge, pair);
 
     if (ivf != ivl) {
       if (orEdge == TopAbs_FORWARD)
-	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier));
       else if (orEdge == TopAbs_REVERSED)
-	structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier));
       else if (orEdge == TopAbs_INTERNAL)
-	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed, nbDomains));
+        structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed));
     }
       
 
   }
-  // Si on a dÝjÒ vu cette Edge et qu`elle n`est pas dÝgÝnÝrÝe on reprend
-  // les points du polygone calculÝs Ò la premiÛre rencontre :
+  // If this Edge has been already checked and it is not degenerated, 
+  // the points of the polygon calculated at the first check are retrieved :
   else
   {
     if (degener)
@@ -1181,14 +997,13 @@
     }
     else
     {
-      // recuperation du polygone:
-      const TColStd_ListOfTransient& L = edges.Find(edge);
-      const Handle(Poly_PolygonOnTriangulation)& P = 
-        *(Handle(Poly_PolygonOnTriangulation)*)&(L.First());
+      // retrieve the polygone:
+      const BRepMesh_PairOfPolygon& pair = edges.Find(edge);
+      const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
       const TColStd_Array1OfInteger& NOD = P->Nodes();
       Handle(TColStd_HArray1OfReal) Par = P->Parameters();
 
-      // creation du nouveau:
+      // creation anew:
       const Standard_Integer nbnodes = NOD.Length();
       TColStd_Array1OfInteger Nodes(1, nbnodes);
       TColStd_Array1OfReal PPar(1, nbnodes);
@@ -1209,7 +1024,7 @@
           {
             const Standard_Real puv = Par->Value(i);
 	    C2d->D0(puv, uv);
-	    v2.Initialize(uv.Coord(), nbDomains, NOD(i), MeshDS_OnCurve);
+	    v2.Initialize(uv.Coord(), NOD(i), MeshDS_OnCurve);
             iv2 = structure->AddNode(v2);
 	
 	    Standard_Integer isv = myvemap.FindIndex(iv2);
@@ -1218,11 +1033,11 @@
 	    PPar(i) = puv;
             
             if (orEdge==TopAbs_FORWARD)
-	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier));
 	    else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier));
             else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed));
 	  
             ivf = iv2;
 	  }
@@ -1256,7 +1071,7 @@
 	    pcos.Perform(P3d, wCurFound);
 	    if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
 	    C2d->D0(wCurFound, uv);
-	    v2.Initialize(uv.Coord(), nbDomains, NOD(i), MeshDS_OnCurve);
+	    v2.Initialize(uv.Coord(), NOD(i), MeshDS_OnCurve);
             iv2 = structure->AddNode(v2);
             
 	    Standard_Integer isv = myvemap.FindIndex(iv2);
@@ -1265,11 +1080,11 @@
 	    PPar(i) = wCurFound;
             
             if (orEdge==TopAbs_FORWARD)
-	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Frontier));
             else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv2, ivf, MeshDS_Frontier));
             else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed, nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivf, iv2, MeshDS_Fixed));
 	  
             ivf = iv2;
 	  }
@@ -1280,11 +1095,11 @@
     
     if (ivf != ivl) {
       if (orEdge == TopAbs_FORWARD) 
-	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier, nbDomains));
+	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Frontier));
       else if (orEdge == TopAbs_REVERSED)
-	structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier, nbDomains));
+	structure->AddLink(BRepMesh_Edge(ivl, ivf, MeshDS_Frontier));
       else if (orEdge == TopAbs_INTERNAL)
-	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed, nbDomains));
+	structure->AddLink(BRepMesh_Edge(ivf, ivl, MeshDS_Fixed));
       }
     }
   }
@@ -1292,17 +1107,17 @@
   P1->Deflection(defedge);
   if (internaledges.IsBound(edge))
   {
-    TColStd_ListOfTransient& L = internaledges.ChangeFind(edge);
+    BRepMesh_PairOfPolygon& pair = internaledges.ChangeFind(edge);
     if (orEdge == TopAbs_REVERSED)
-      L.Append(P1);
+      pair.Append(P1);
     else
-      L.Prepend(P1);
+      pair.Prepend(P1);
   }
   else
   {
-    TColStd_ListOfTransient L1;
-    L1.Append(P1);
-    internaledges.Bind(edge, L1);
+    BRepMesh_PairOfPolygon pair1;
+    pair1.Append(P1);
+    internaledges.Bind(edge, pair1);
   }
 }
 
@@ -1369,7 +1184,7 @@
 
 	if (sameUV)
         {
-	  // 1. est-ce vraiment sameUV sans etre denegere
+	  // 1. is it really sameUV without being degenerated
 	  gp_Pnt2d uvF, uvL;
 	  C2d->D0(first, uvF);
 	  C2d->D0(last, uvL);
@@ -1417,7 +1232,7 @@
 	} 
 	NewNodInStruct(1) = ipf;
 	theUV = FindUV(pBegin, uvFirst, ipf, gFace, mindist);
-	BRepMesh_Vertex vf(theUV,nbDomains,ipf,MeshDS_Frontier);
+	BRepMesh_Vertex vf(theUV,ipf,MeshDS_Frontier);
         iv1 = structure->AddNode(vf);
         isv1 = myvemap.FindIndex(iv1);
         if (isv1 == 0) isv1 = myvemap.Add(iv1);
@@ -1454,7 +1269,7 @@
         }
 	NewNodInStruct(nbnodes) = ipl;
 	theUV = FindUV(pEnd, uvLast, ipl, gFace, mindist);
-	BRepMesh_Vertex vl(theUV,nbDomains,ipl,MeshDS_Frontier);
+	BRepMesh_Vertex vl(theUV,ipl,MeshDS_Frontier);
         
         ivl = structure->AddNode(vl);
         isvl = myvemap.FindIndex(ivl);
@@ -1478,7 +1293,7 @@
 	    NewNodInStruct(i) = nbLocat;
             // Record 2d point
 	    uv = C2d->Value(Param->Value(i));
-	    v.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_Frontier);
+	    v.Initialize(uv.Coord(), nbLocat, MeshDS_Frontier);
             iv2 = structure->AddNode(v);
             isv = myvemap.FindIndex(iv2);
             if (isv == 0) isv = myvemap.Add(iv2);
@@ -1486,22 +1301,22 @@
             
             //add links
             if (orEdge == TopAbs_FORWARD)
-              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier,nbDomains));
+              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed));
 	    iv1 = iv2;  
 	  }
           
           // last point
           if (iv1 != ivl) {
 	    if (orEdge == TopAbs_FORWARD)
-	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed));
 	  }
 
           
@@ -1541,7 +1356,7 @@
 	    pcos.Perform(P3d, wCurFound);
 	    if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
 	    C2d->D0(wCurFound, uv);
-	    v.Initialize(uv.Coord(), nbDomains, nbLocat, MeshDS_Frontier);
+	    v.Initialize(uv.Coord(), nbLocat, MeshDS_Frontier);
             iv2 = structure->AddNode(v);
             isv = myvemap.FindIndex(iv2);
             if (isv == 0) isv = myvemap.Add(iv2);
@@ -1550,22 +1365,22 @@
             
             //add links
             if (orEdge == TopAbs_FORWARD)
-              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier,nbDomains));
+              structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv2,iv1,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,iv2,MeshDS_Fixed));
 	    iv1 = iv2;              
 	  }
           
           // last point
           if (iv1 != ivl) {
 	    if (orEdge == TopAbs_FORWARD)
-	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_REVERSED)
-	      structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(ivl,iv1,MeshDS_Frontier));
 	    else if (orEdge == TopAbs_INTERNAL)
-	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed,nbDomains));
+	      structure->AddLink(BRepMesh_Edge(iv1,ivl,MeshDS_Fixed));
 	  }
 	}
 
@@ -1574,25 +1389,25 @@
 	P1->Deflection(defedge);
 	if (internaledges.IsBound(edge))
         {
-	  TColStd_ListOfTransient& L = internaledges.ChangeFind(edge);
+	  BRepMesh_PairOfPolygon& pair = internaledges.ChangeFind(edge);
           if (edge.Orientation() == TopAbs_REVERSED)
-            L.Append(P1);
+            pair.Append(P1);
           else
-            L.Prepend(P1);
+            pair.Prepend(P1);
 	}
 	else
         {
-	  TColStd_ListOfTransient L1;
-	  L1.Append(P1);
-	  internaledges.Bind(edge, L1);
+	  BRepMesh_PairOfPolygon pair1;
+	  pair1.Append(P1);
+	  internaledges.Bind(edge, pair1);
 	}
 
 	Handle(Poly_PolygonOnTriangulation) P2 =
           new Poly_PolygonOnTriangulation(NewNodInStruct, Param->Array1());
 	P2->Deflection(defedge);
-	TColStd_ListOfTransient L;
-	L.Append(P2);
-	edges.Bind(edge, L);
+	BRepMesh_PairOfPolygon pair;
+	pair.Append(P2);
+	edges.Bind(edge, pair);
 
 	found = Standard_True;
       }
@@ -1615,1383 +1430,150 @@
   return found;
 }
 
+//=======================================================================
+//function : output
+//purpose  : 
+//=======================================================================
+Standard_Integer BRepMesh_FastDiscret::NbTriangles() const
+{
+  return structure->NbElements();
+}
+
+//=======================================================================
+//function : Triangle
+//purpose  : 
+//=======================================================================
+
+const BRepMesh_Triangle& BRepMesh_FastDiscret::Triangle
+  (const Standard_Integer Index) const
+{
+  return structure->GetElement(Index);
+}
+
+//=======================================================================
+//function : NbEdges
+//purpose  : 
+//=======================================================================
+
+Standard_Integer BRepMesh_FastDiscret::NbEdges() const
+{
+  return structure->NbLinks();
+}
+
+//=======================================================================
+//function : Edge
+//purpose  : 
+//=======================================================================
+
+const BRepMesh_Edge& BRepMesh_FastDiscret::Edge(const Standard_Integer Index) const
+{
+  return structure->GetLink(Index);
+}
+
+//=======================================================================
+//function : NbVertices
+//purpose  : 
+//=======================================================================
+
+Standard_Integer BRepMesh_FastDiscret::NbVertices() const
+{
+  return structure->NbNodes();
+}
 
+//=======================================================================
+//function : Vertex
+//purpose  : 
+//=======================================================================
+
+const BRepMesh_Vertex& BRepMesh_FastDiscret::Vertex
+  (const Standard_Integer Index) const
+{
+  return structure->GetNode(Index);
+}
 
 //=======================================================================
-//function : InternalVertices
+//function : Pnt
 //purpose  : 
 //=======================================================================
-void BRepMesh_FastDiscret::InternalVertices ( const Handle(BRepAdaptor_HSurface)& caro,
-                                          BRepMesh_ListOfVertex&              InternalV,
-                                          const Standard_Real                 defface,
-                                          const BRepMesh_Classifier&          classifier)
+
+const gp_Pnt& BRepMesh_FastDiscret::Pnt(const Standard_Integer Index) const
 {
-  BRepMesh_Vertex newV;
-  gp_Pnt2d p2d;
-  gp_Pnt p3d;
-  MeshShape_ListOfSurfacePoint pntsOnSurf;
+  return Location3d(structure->GetNode(Index).Location3d());
+}
+
+//=======================================================================
+//function : VerticesOfDomain
+//purpose  : 
+//=======================================================================
+
+void BRepMesh_FastDiscret::VerticesOfDomain(MeshDS_MapOfInteger&  Indices) const 
+{ 
+  Indices.Clear();
+  
+  // recuperate from the map of edges.
+  const MeshDS_MapOfInteger& edmap = structure->LinkOfDomain();
+
+  // iterator on edges.
+  MeshDS_MapOfInteger::Iterator iter(edmap);
   
-  // travail suivant le type de surface
+  Standard_Integer ind_edge;
+  for (iter.Reset(); iter.More(); iter.Next()) {
+    ind_edge = iter.Key();
+    const BRepMesh_Edge& Ed = Edge(ind_edge);
+    Indices.Add(Ed.FirstNode());
+    Indices.Add(Ed.LastNode());
+  }
+}
 
-  const BRepAdaptor_Surface& BS = *(BRepAdaptor_Surface*)&(caro->Surface());
-  GeomAbs_SurfaceType thetype = caro->GetType();
+//=======================================================================
+//function : FindUV
+//purpose  : 
+//=======================================================================
 
-  if (thetype == GeomAbs_Plane && !classifier.NaturalRestriction())
+gp_XY BRepMesh_FastDiscret::FindUV(const TopoDS_Vertex&   V,
+			       const gp_Pnt2d&        XY, 
+			       const Standard_Integer ip,
+			       const Handle(BRepAdaptor_HSurface)& S,
+			       const Standard_Real mindist) 
+{
+  gp_XY theUV;
+  if (mylocation2d.IsBound(ip))
   {
-    // rajout d`un seul point au milieu.
-    const Standard_Real U = 0.5*(myumin+myumax);
-    const Standard_Real V = 0.5*(myvmin+myvmax);
-    if (classifier.Perform(gp_Pnt2d(U, V)) == TopAbs_IN)
-    {
-      // Record 3d point
-      BRepMesh_GeomTool::D0(caro, U, V, p3d);
-      nbLocat++;
-      Location3d.Bind(nbLocat, p3d);
-      // Record 2d point
-      p2d.SetCoord((U-myumin)/deltaX, (V-myvmin)/deltaY);
-      newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-      InternalV.Append(newV);
-    }
-  }
-  else if (thetype == GeomAbs_Sphere)
-  {
-    gp_Sphere S = BS.Sphere();
-    const Standard_Real R = S.Radius();
-
-    // Calculate parameters for iteration in V direction
-    Standard_Real Dv = 1.0 - (defface/R);
-    if (Dv < 0.0) Dv = 0.0;
-    Standard_Real oldDv = 2.0 * ACos (Dv);
-    Dv  =  .7 * oldDv; //.7 ~= sqrt(2.) - Dv is hypotenuse of triangle when oldDv is legs
-    const Standard_Real sv = myvmax - myvmin;
-    Dv = sv/((Standard_Integer)(sv/Dv) + 1);
-    const Standard_Real pasvmax = myvmax-Dv*0.5;
-
-    //Du can be defined from relation: 2*r*Sin(Du/2) = 2*R*Sin(Dv/2), r = R*Cos(v)
-    //here approximate relation r*Du = R*Dv is used
-    
-    Standard_Real Du, pasu, pasv; //, ru;
-    const Standard_Real su = myumax-myumin;
-    Standard_Boolean Shift = Standard_False;
-    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv)
+    BRepMesh_ListOfXY& L = mylocation2d.ChangeFind(ip);
+    theUV = L.First();
+    if (L.Extent() != 1)
     {
-      // Calculate parameters for iteration in U direction
-      // 1.-.365*pasv*pasv is simple approximation of Cos(pasv)
-      // with condition that it gives ~.1 when pasv = pi/2
-      Du = Dv/(1.-.365*pasv*pasv);
-      Du = su/((Standard_Integer)(su/Du) + 1);
-      Shift = !Shift;
-      const Standard_Real d = (Shift)? Du*.5 : 0.;
-      const Standard_Real pasumax = myumax-Du*0.5 + d;
-      for (pasu = myumin + Du - d; pasu < pasumax; pasu += Du)
+      BRepMesh_ListIteratorOfListOfXY it(L);
+      it.Next();
+      Standard_Real dd, dmin = XY.Distance(gp_Pnt2d(theUV));
+      for (; it.More(); it.Next())
       {
-	if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
+	dd = XY.Distance(gp_Pnt2d(it.Value()));
+	if (dd < dmin)
         {
-          // Record 3d point
-#ifdef DEB_MESH_CHRONO
-	  D0Internal++;
-#endif
-	  ElSLib::D0(pasu, pasv, S, p3d);
-	  nbLocat++;
-	  Location3d.Bind(nbLocat, p3d);
-          // Record 2d point
-	  p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
-	  newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-	  InternalV.Append(newV);
+	  theUV = it.Value();
+	  dmin = dd;
 	}
       }
     }
-  }
-  else if (thetype == GeomAbs_Cylinder)
-  {
-    gp_Cylinder S = BS.Cylinder();
-    const Standard_Real R = S.Radius();
 
-    // Calculate parameters for iteration in U direction
-    Standard_Real Du = 1.0 - (defface/R);
-    if (Du < 0.0) Du = 0.0;
-    Du = 2.0 * ACos (Du);
-    if (Du > angle) Du = angle;
-    const Standard_Real su = myumax - myumin;
-    const Standard_Integer nbU = (Standard_Integer)(su/Du);
-    Du = su/(nbU+1);
-
-    // Calculate parameters for iteration in V direction
-    const Standard_Real sv = myvmax - myvmin;
-    Standard_Integer nbV = (Standard_Integer)( nbU*sv/(su*R) );
-    nbV = Min(nbV, 100*nbU);
-    Standard_Real Dv = sv/(nbV+1);
-
-    Standard_Real pasu, pasv, pasvmax = myvmax-Dv*0.5, pasumax = myumax-Du*0.5;
-    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv) {
-      for (pasu = myumin + Du; pasu < pasumax; pasu += Du) {
-	if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
-        {
-          // Record 3d point
-	  ElSLib::D0(pasu, pasv, S, p3d);
-          nbLocat++;
-          Location3d.Bind(nbLocat, p3d);
-          // Record 2d point
-	  p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
-	  newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-          InternalV.Append(newV);
-	}
-      }
+    const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(V), mindist);
+
+    const Standard_Real Utol2d = .5 * (S->LastUParameter() - S->FirstUParameter());
+    const Standard_Real Vtol2d = .5 * (S->LastVParameter() - S->FirstVParameter());
+
+    const gp_Pnt p1 = S->Value(theUV.X(), theUV.Y());
+    const gp_Pnt p2 = S->Value(XY.X(), XY.Y());
+
+    if (Abs(theUV.X() - XY.X()) > Utol2d ||
+        Abs(theUV.Y() - XY.Y()) > Vtol2d ||
+	!p1.IsEqual(p2, tol))
+    {
+      theUV = XY.Coord();
+      L.Append(theUV);
     }
   }
-  else if (thetype == GeomAbs_Cone)
-  {
-    Standard_Real R, RefR, SAng;
-    gp_Cone C = BS.Cone();
-    RefR = C.RefRadius();
-    SAng = C.SemiAngle();
-    R = Max(Abs(RefR+myvmin*Sin(SAng)), Abs(RefR+myvmax*Sin(SAng)));
-    Standard_Real Du, Dv, pasu, pasv;
-    Du = Max(1.0e0 - (defface/R),0.0e0);
-    Du  = (2.0 * ACos (Du));
-    Standard_Integer nbU = (Standard_Integer) ( (myumax-myumin)/Du );
-    Standard_Integer nbV = (Standard_Integer) ( nbU*(myvmax-myvmin)/((myumax-myumin)*R) );
-    Du = (myumax-myumin)/(nbU+1);
-    Dv = (myvmax-myvmin)/(nbV+1);
-
-    Standard_Real pasvmax = myvmax-Dv*0.5, pasumax = myumax-Du*0.5;
-    for (pasv = myvmin + Dv; pasv < pasvmax; pasv += Dv) {
-      for (pasu = myumin + Du; pasu < pasumax; pasu += Du) {
-	if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
-        {
-          // Record 3d point
-	  ElSLib::D0(pasu, pasv, C, p3d);
-          nbLocat++;
-          Location3d.Bind(nbLocat, p3d);
-          // Record 2d point
-          p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
-	  newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-          InternalV.Append(newV);
-	}
-      }
-    }
-  }
-  else if (thetype == GeomAbs_Torus)
-  {
-    gp_Torus T = BS.Torus();
-
-    Standard_Boolean insert;
-    Standard_Integer i, j, ParamULength, ParamVLength;
-    Standard_Real pp, pasu, pasv;
-    Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
-
-    TColStd_SequenceOfReal ParamU, ParamV;
-
-    Standard_Real Du, Dv;//, pasu, pasv;
-    Dv = Max(1.0e0 - (defface/r),0.0e0) ;
-    Standard_Real oldDv = 2.0 * ACos (Dv);
-    oldDv = Min(oldDv, angle);
-    Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
-    Dv = oldDv;
-    
-    Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
-    Dv = (myvmax-myvmin)/(nbV+1);
-    Standard_Real ru = R + r;
-    if (ru > 1.e-16)
-    {
-      Du  = 2.0 * ACos(Max(1.0 - (defface/ru),0.0));
-      if (angle < Du) Du = angle;
-      Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
-      if(aa < gp::Resolution())
-        return; 
-      Du *= Min(oldDv, Du) / aa;
-    }
-    else Du = Dv;     
-    
-    Standard_Integer nbU = Max((Standard_Integer)((myumax-myumin)/Du), 2);
-    nbU = Max(nbU , (int)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
-    Du = (myumax-myumin)/(nbU+1);
-    
-    if (R < r)
-    {
-      // comme on recupere les points des edges.
-      // dans ce cas, les points ne sont pas representatifs.
-            
-      //-- On choisit DeltaX et DeltaY de facon a ce qu on ne saute pas 
-      //-- de points sur la grille
-      for (i = 0; i <= nbU; i++) ParamU.Append(myumin + i* Du);
-    }//R<r
-    else //U if R > r
-    {
-      //--ofv: U
-      // Number of mapped U parameters
-      const Standard_Integer LenU = myUParam.Extent();
-      // Fill array of U parameters
-      TColStd_Array1OfReal Up(1,LenU);
-      for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
-      
-      // Calculate DU, sort array of parameters
-      Standard_Real aDU = FUN_CalcAverageDUV(Up,LenU);
-      aDU = Max(aDU, Abs(myumax -  myumin) / (Standard_Real) nbU / 2.);
-      Standard_Real dUstd = Abs(myumax -  myumin) / (Standard_Real) LenU;
-      if (aDU > dUstd) dUstd = aDU;
-      // Add U parameters
-      for (j = 1; j <= LenU; j++)
-      {
-        pp = Up(j);
-        insert = Standard_True;
-        ParamULength = ParamU.Length();
-        for (i = 1; i <= ParamULength && insert; i++)
-        {
-          insert = (Abs(ParamU.Value(i)-pp) > (0.5*dUstd));
-        }
-        if (insert) ParamU.Append(pp);
-      }
-    } 
-
-    //--ofv: V
-    // Number of mapped V parameters
-    const Standard_Integer LenV = myVParam.Extent();
-    // Fill array of V parameters
-    TColStd_Array1OfReal Vp(1,LenV);
-    for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
-    // Calculate DV, sort array of parameters
-    Standard_Real aDV = FUN_CalcAverageDUV(Vp,LenV);
-    aDV = Max(aDV, Abs(myvmax -  myvmin) / (Standard_Real) nbV / 2.);
-
-    Standard_Real dVstd = Abs(myvmax -  myvmin) / (Standard_Real) LenV;
-    if (aDV > dVstd) dVstd = aDV;
-    // Add V parameters
-    for (j = 1; j <= LenV; j++)
-    {
-      pp = Vp(j);
-
-      insert = Standard_True;
-      ParamVLength = ParamV.Length();
-      for (i = 1; i <= ParamVLength && insert; i++)
-      {
-        insert = (Abs(ParamV.Value(i)-pp) > (dVstd*2./3.));
-      }
-      if (insert) ParamV.Append(pp);
-    }
-
-    Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
-    Standard_Real umin = myumin+deltaY*0.1;
-    Standard_Real vmin = myvmin+deltaX*0.1;
-    Standard_Real umax = myumax-deltaY*0.1;
-    Standard_Real vmax = myvmax-deltaX*0.1;
-
-    for (i = 1; i <= Lu; i++)
-    {
-      pasu = ParamU.Value(i);
-      if (pasu >= umin && pasu < umax)
-      {
-	for (j = 1; j <= Lv; j++)
-        {
-	  pasv = ParamV.Value(j);
-	  if (pasv >= vmin && pasv < vmax)
-          {
-	    if (classifier.Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
-            {
-              // Record 3d point
-	      ElSLib::D0(pasu, pasv, T, p3d);
-	      nbLocat++;
-	      Location3d.Bind(nbLocat, p3d);
-              // Record 2d point
-	      p2d.SetCoord((pasu-myumin)/deltaX, (pasv-myvmin)/deltaY);
-	      newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-	      InternalV.Append(newV);
-	    }
-	  }
-	}
-      }
-    }
-  }
-  else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
-  {
-    Standard_Integer i, j;
-    Standard_Real pp;
-
-    // Sort sequence of U parameters
-    TColStd_SequenceOfReal ParamU;
-    ParamU.Append(myumin); ParamU.Append(myumax);
-    Standard_Integer ParamULength = 2;
-    const Standard_Integer LenU = myUParam.Extent();
-    const Standard_Real ddu = 0.02*(myumax-myumin);
-    
-    Standard_Real U1, U2;
-
-    for (j = 1; j <= LenU; j++)
-    {
-      pp = myUParam(j);
-
-      U1 = ParamU.Value(1);
-      for (i = 2; i <= ParamULength; i++)
-      {
-	U2 = ParamU.Value(i);
-	if (pp < U2)
-        {
-	  if ((U2-pp) < ddu)
-          {
-	    ParamU(i) = pp;
-	  }
-	  else if ((pp-U1) > ddu)
-          {
-	    ParamU.InsertBefore(i, pp);
-	    ParamULength++;
-	  }
-	  else if (i != 2)
-          {
-	    ParamU(i-1) = pp;
-	  }
-	  break;
-	}
-        U1 = U2;
-      }
-    }
-    
-    // Sort sequence of V parameters
-    TColStd_SequenceOfReal ParamV;
-    ParamV.Append(myvmin); ParamV.Append(myvmax);
-    Standard_Integer ParamVLength = 2;
-    const Standard_Integer LenV = myVParam.Extent();
-    const Standard_Real ddv = 0.02*(myvmax-myvmin);
-
-    Standard_Real V1, V2;
-    for (j = 1; j <= LenV; j++)
-    {
-      pp = myVParam(j);
-
-      V1 = ParamV.Value(1);
-      for (i = 2; i <= ParamVLength; i++)
-      {
-	V2 = ParamV.Value(i);
-	if (pp < V2)
-        {
-	  if ((V2-pp) < ddv)
-          {
-	    ParamV(i) = pp;
-	  }
-	  else if ((pp-V1) > ddv)
-          {
-	    ParamV.InsertBefore(i, pp);
-	    ParamVLength++;
-	  }
-	  else if (i != 2)
-          {
-	    ParamV(i-1) = pp;
-	  }
-	  break;
-	}
-        V1 = V2;
-      }
-    }
-
-    // controle des isos U et insertion eventuelle:
-
-    Handle(Geom_Surface) B;
-    if (thetype == GeomAbs_BezierSurface) {
-      B = BS.Bezier();
-    }
-    else {
-      B = BS.BSpline();
-    }
-
-    gp_Pnt P1, P2, PControl;
-    Standard_Real u, v, dist;
-
-    // Insert V parameters by deflection criterion
-    for (i = 1; i <= ParamULength; i++)
-    {
-      Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
-      V1 = ParamV.Value(1);
-      P1 = IsoU->Value(V1);
-      for (j = 2; j <= ParamVLength;)
-      {
-	V2 = ParamV.Value(j);
-	P2 = IsoU->Value(V2);
-	v = 0.5*(V1+V2);
-	PControl = IsoU->Value(v);
-	if (P1.SquareDistance(P2) > 1.e-5)
-        {
-	  gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
-	  dist = L.Distance(PControl);
-	}
-	else
-        {
-	  dist = P1.Distance(PControl);
-	}
-        if (dist > defface)
-        {
-          // insertion 
-          ParamV.InsertBefore(j, v);
-          ParamVLength++;
-        }
-        else
-        {
-          V1 = V2;
-          P1 = P2;
-          j++;
-        }
-      }
-    }
-
-    for (i = 2; i < ParamVLength; i++)
-    {
-      v = ParamV.Value(i);
-      Handle(Geom_Curve) IsoV = B->VIso(v);
-      U1 = ParamU.Value(1);
-      P1 = IsoV->Value(U1);
-      for (j = 2; j <= ParamULength;)
-      {
-	U2 = ParamU.Value(j);
-	P2 = IsoV->Value(U2);
-	u = 0.5*(U1+U2);
-	PControl = IsoV->Value(u);
-	if (P1.SquareDistance(P2) > 1.e-5)
-        {
-	  gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
-	  dist = L.Distance(PControl);
-	}
-	else
-        {
-	  dist = P1.Distance(PControl);
-	}
-        if (dist > defface)
-        {
-          // insertion 
-          ParamU.InsertBefore(j, u);
-          ParamULength++;
-        }
-        else
-        {
-          U1 = U2;
-          P1 = P2;
-          if (j < ParamULength)
-          {
-            // Classify intersection point
-            if (classifier.Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
-            {
-              // Record 3d point
-              nbLocat++;
-              Location3d.Bind(nbLocat, P1);
-              // Record 2d point
-              p2d.SetCoord((U1-myumin)/deltaX, (v-myvmin)/deltaY);
-              newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-              InternalV.Append(newV);
-            }
-          }
-          j++;
-        }
-      }
-    }
-  }
-  else
-  {
-    const Standard_Real theangle = 0.35;
-
-    Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
-    Adaptor3d_IsoCurve tabu[10], tabv[10];
-
-    TColStd_SequenceOfReal ParamU, ParamV;
-    Standard_Real u, v, du, dv;
-    Standard_Integer iu, iv;
-    Standard_Real f, l;
-
-    du = (myumax-myumin) / (nbpointsU+1); dv = (myvmax-myvmin) / (nbpointsV+1);
-
-    for (iu = 1; iu <= nbpointsU; iu++) {
-      u = myumin + iu*du;
-      tabu[iu-1].Load(caro);
-      tabu[iu-1].Load(GeomAbs_IsoU, u);
-    }
-
-    for (iv = 1; iv <= nbpointsV; iv++) {
-      v = myvmin + iv*dv;
-      tabv[iv-1].Load(caro);
-      tabv[iv-1].Load(GeomAbs_IsoV, v);
-    }
-
-    Standard_Integer imax = 1, MaxV = 0;
-    
-    GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU];
-
-    for (i = 0; i <= nbpointsU-1; i++) {
-      f = Max(myvmin, tabu[i].FirstParameter());
-      l = Min(myvmax, tabu[i].LastParameter());
-      GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, theangle, 0.7*defface, 2);
-      tabGU[i] = theDeflection;
-      if (tabGU[i].NbPoints() > MaxV) {
-	MaxV = tabGU[i].NbPoints();
-	imax = i;
-      }
-    }
-    
-    // recuperation du tableau de parametres V:
-    Standard_Integer NV = tabGU[imax].NbPoints();
-    for (i = 1; i <= NV; i++) {
-      ParamV.Append(tabGU[imax].Parameter(i));
-    }
-    delete [] tabGU;
-
-    imax = 1;
-    Standard_Integer MaxU = 0;
-
-    GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV];
-
-    for (i = 0; i <= nbpointsV-1; i++) {
-      f = Max(myumin, tabv[i].FirstParameter());
-      l = Min(myumax, tabv[i].LastParameter());
-      GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, theangle, 0.7*defface, 2);
-      tabGV[i] = thedeflection2;
-      if (tabGV[i].NbPoints() > MaxU) {
-	MaxU = tabGV[i].NbPoints();
-	imax = i;
-      }
-    }
-    
-    // recuperation du tableau de parametres U:
-    Standard_Integer NU = tabGV[imax].NbPoints();
-    for (i = 1; i <= NU; i++) {
-      ParamU.Append(tabGV[imax].Parameter(i));
-    }
-    delete [] tabGV;
-    
-    if (ParamU.Length() == 2) {
-      ParamU.InsertAfter(1, (myumax+myumin)*0.5);
-    }
-    if (ParamV.Length() == 2) {
-      ParamV.InsertAfter(1, (myvmax+myvmin)*0.5);
-    }
-    
-    TColStd_SequenceOfReal InsertV, InsertU;
-    gp_Pnt P1;
-
-    Adaptor3d_IsoCurve IsoV;
-    IsoV.Load(caro);
-
-    Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
-    
-    for (i = 2; i < Lv; i++) {
-      v = ParamV.Value(i);
-      IsoV.Load(GeomAbs_IsoV, v);
-      for (j = 2; j < Lu; j++) {
-	u = ParamU.Value(j);
-	if (classifier.Perform(gp_Pnt2d(u, v)) == TopAbs_IN)
-        {
-          // Record 3d point
-	  P1 = IsoV.Value(u);
-	  nbLocat++;
-          Location3d.Bind(nbLocat, P1);
-          // Record 2d point
-	  p2d.SetCoord((u-myumin)/deltaX, (v-myvmin)/deltaY);
-	  newV.Initialize(p2d.XY(), nbDomains, nbLocat, MeshDS_Free);
-	  InternalV.Append(newV); 
-	}
-      }
-    } 
-  }
-#ifdef DEB_MESH_CHRONO
-  chInternal.Stop();
-#endif
-}
-
-//=======================================================================
-//function : InternalVerticesForRectangle
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscret::
-  InternalVerticesForRectangle(const Handle(BRepAdaptor_HSurface)& caro,
-			       BRepMesh_ListOfVertex&              InternalV,
-			       const TColStd_ListOfReal&           U1Params,
-			       const TColStd_ListOfReal&           U2Params,
-			       const TColStd_ListOfReal&           V1Params,
-			       const TColStd_ListOfReal&           V2Params)
-{
-  BRepMesh_Vertex newV;
-  gp_XY p2d;
-  gp_Pnt p3d;
-
-
-  Standard_Integer i, j, nbu = U1Params.Extent(), 
-                         nbv = V1Params.Extent();
-  TColStd_Array1OfReal u1prs(1, nbu), u2prs(1, nbu), 
-                       v1prs(1, nbv), v2prs(1,nbv);
-  TColStd_ListIteratorOfListOfReal Itu1(U1Params), Itu2(U2Params), 
-                                   Itv1(V1Params), Itv2(V2Params);
-  for(i = 1; Itu1.More(); Itu1.Next(), Itu2.Next(),++i) {
-    u1prs(i) = Itu1.Value();
-    u2prs(i) = Itu2.Value();
-  }
-  for(i = 1; Itv1.More(); Itv1.Next(), Itv2.Next(),++i) {
-    v1prs(i) = Itv1.Value();
-    v2prs(i) = Itv2.Value();
-  }
-  
-  SortTools_ShellSortOfReal aSorter;
-  TCollection_CompareOfReal aCompare;
-  aSorter.Sort(u1prs, aCompare);
-  aSorter.Sort(u2prs, aCompare);
-  aSorter.Sort(v1prs, aCompare);
-  aSorter.Sort(v2prs, aCompare);
-
-  gp_XY f00(u1prs(1), v1prs(1));
-  gp_XY f01(u2prs(1), v1prs(nbv));
-  gp_XY f10(u1prs(nbu), v2prs(1));
-  gp_XY f11(u2prs(nbu), v2prs(nbv));
-
-  Standard_Real cu = 1./(nbu-1), cv = 1./(nbv-1);
-
-  //Linear Coons-Gordon interpolation for inner points
-  Standard_Real x, y, x1, y1;
-  for(i = 2, x = cu; i < nbu; ++i, x += cu) {
-    x1 = 1.-x;
-    gp_XY fx0(u1prs(i), v1prs(1));
-    gp_XY fx1(u2prs(i), v1prs(nbv));
- 
-    for(j = 2, y = cv; j < nbv; ++j, y += cv) {
-      y1 = 1.-y;
-
-      gp_XY f0y(u1prs(1), v1prs(j));
-      gp_XY f1y(u1prs(nbu), v2prs(j));
-		
-      gp_XY fxy = y1*fx0 + y*fx1 + x1*f0y + x*f1y -
-                  (x1*(y1*f00 + y*f01) + x*(y1*f10 + y*f11));
-
-      // Record 3d point
-      BRepMesh_GeomTool::D0(caro, fxy.X(), fxy.Y(), p3d);
-      nbLocat++;
-      Location3d.Bind(nbLocat, p3d);
-      // Record 2d point
-      p2d.SetCoord((fxy.X()-myumin)/deltaX, (fxy.Y()-myvmin)/deltaY);
-      newV.Initialize(p2d, nbDomains, nbLocat, MeshDS_Free);
-      InternalV.Append(newV);
-    }
-  }
-      
-
-
-}
-
-
-//=======================================================================
-//function : Control
-//purpose  : 
-//=======================================================================
-Standard_Real BRepMesh_FastDiscret::Control
-  (const Handle(BRepAdaptor_HSurface)& caro,
-   const Standard_Real                 defface,
-   BRepMesh_ListOfVertex&              InternalV,
-   TColStd_ListOfInteger&              badTriangles,
-   TColStd_ListOfInteger&              nulTriangles,
-   BRepMesh_Delaun&                    trigu,
-   const Standard_Boolean              isfirst)
-{
-  //IMPORTANT: Constants used in calculations
-  const Standard_Real MinimalArea2d = 1.e-9;
-  const Standard_Real MinimalSqLength3d = 1.e-12;
-
-  // Define the number of iterations
-  Standard_Integer myNbIterations = 11;
-  const Standard_Integer nbPasses = (isfirst? 1 : myNbIterations);
-
-  // Initialize stop condition
-  Standard_Boolean allDegenerated = Standard_False;
-  Standard_Integer nbInserted = 1;
-
-  // Create map of links to skip already processed
-  Standard_Integer nbtriangles;
-
-  nbtriangles = structure->ElemOfDomain(nbDomains).Extent();
-  if (nbtriangles <= 0) return -1.0;
-  MeshShape_MapOfCouple theCouples(3*nbtriangles);
-
-  gp_XY mi2d;
-  gp_XYZ vecEd1, vecEd2, vecEd3;
-  gp_Pnt pDef;
-  Standard_Real dv, defl, maxdef;
-  Standard_Integer pass, nf, nl;
-  Standard_Integer v1, v2, v3, e1, e2, e3;
-  Standard_Boolean o1, o2, o3;
-  Standard_Boolean m1, m2, m3;
-  BRepMesh_Vertex InsVertex;
-  Standard_Boolean caninsert;
-
-  // Perform refinement passes
-  for (pass = 1; pass <= nbPasses && nbInserted && !allDegenerated; pass++)
-  {
-    InternalV.Clear();
-    badTriangles.Clear();
-    
-    // Reset stop condition
-    allDegenerated = Standard_True;
-    nbInserted = 0;
-    maxdef = -1.0;
-
-    // Do not insert nodes in last pass in non-SharedMode
-    caninsert = (WithShare || pass < nbPasses);
-
-    // Read mesh size
-    nbtriangles = structure->ElemOfDomain(nbDomains).Extent();
-    if (nbtriangles <= 0) break;
-
-    // Iterate on current triangles
-    BRepMesh_ListOfVertex InsVertices;
-    TColStd_MapOfInteger usedEdges;
-   // TColStd_MapIteratorOfMapOfInteger triDom;
-    MeshDS_MapOfInteger::Iterator triDom;
-    const MeshDS_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
-    triDom.Initialize(TriMap);
-    Standard_Integer aNbPnt = 0;
-    for (; triDom.More(); triDom.Next())
-    {
-      const BRepMesh_Triangle& curTri=Triangle(triDom.Key());
-      if (curTri.Movability()==MeshDS_Deleted) continue;
-      curTri.Edges(e1, e2, e3, o1, o2, o3);
-      
-      const BRepMesh_Edge& edg1=Edge(e1);
-      const BRepMesh_Edge& edg2=Edge(e2);
-      const BRepMesh_Edge& edg3=Edge(e3);
-      
-      m1 = (edg1.Movability() == MeshDS_Frontier);
-      m2 = (edg2.Movability() == MeshDS_Frontier);
-      m3 = (edg3.Movability() == MeshDS_Frontier);
-      if (o1) {
-	v1=edg1.FirstNode();
-	v2=edg1.LastNode();
-      }
-      else {
-	v1=edg1.LastNode();
-	v2=edg1.FirstNode();
-      }
-      if (o2)
-	v3=edg2.LastNode();
-      else
-	v3=edg2.FirstNode();
-
-      const BRepMesh_Vertex& vert1=Vertex(v1);
-      const BRepMesh_Vertex& vert2=Vertex(v2);
-      const BRepMesh_Vertex& vert3=Vertex(v3);
-
-      const gp_XYZ& p1=Location3d(vert1.Location3d()).Coord();
-      const gp_XYZ& p2=Location3d(vert2.Location3d()).Coord();
-      const gp_XYZ& p3=Location3d(vert3.Location3d()).Coord();
-
-      vecEd1 = p2 - p1;
-      vecEd2 = p3 - p2;
-      vecEd3 = p1 - p3;
-
-      // Check for degenerated triangle
-      if (vecEd1.SquareModulus() < MinimalSqLength3d ||
-          vecEd2.SquareModulus() < MinimalSqLength3d ||
-          vecEd3.SquareModulus() < MinimalSqLength3d) 
-      {
-        nulTriangles.Append(triDom.Key());
-        continue;
-      }
-
-      allDegenerated = Standard_False;
-
-      gp_XY xy1(vert1.Coord().X()*deltaX+myumin,vert1.Coord().Y()*deltaY+myvmin);
-      gp_XY xy2(vert2.Coord().X()*deltaX+myumin,vert2.Coord().Y()*deltaY+myvmin);
-      gp_XY xy3(vert3.Coord().X()*deltaX+myumin,vert3.Coord().Y()*deltaY+myvmin);
-
-      // Check triangle area in 2d
-      if (Abs((xy2-xy1)^(xy3-xy1)) < MinimalArea2d)
-      {
-       	nulTriangles.Append(triDom.Key());
-        continue;
-      }
-
-      // Check triangle normal
-      gp_XYZ normal(vecEd1^vecEd2);
-      dv = normal.Modulus();
-      if (dv < Precision::Confusion())
-      {
-       	nulTriangles.Append(triDom.Key());
-        continue;
-      }
-      normal /= dv;
-
-      // Check deflection on triangle
-      mi2d = (xy1+xy2+xy3)/3.0;
-      caro->D0(mi2d.X(), mi2d.Y(), pDef);
-      defl = Abs(normal*(pDef.XYZ()-p1));
-      if (defl > maxdef) maxdef = defl;
-      if (defl > defface)
-      {
-        if (isfirst) break;
-        if (caninsert)
-        {
-          // Record new vertex
-          aNbPnt++;
-          nbLocat++;
-          Location3d.Bind(nbLocat,pDef);
-          mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
-          InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
-          InternalV.Append(InsVertex);
-        }
-        badTriangles.Append(triDom.Key());
-      }
-      
-      if (!m2) // Not a boundary
-      {
-        // Check if this link was already processed
-        if (v2 < v3) { nf = v2; nl = v3; } else { nf = v3; nl = v2; }
-        if (theCouples.Add(MeshShape_Couple(nf,nl)))
-        {
-          // Check deflection on edge 1
-          mi2d = (xy2+xy3)*0.5;
-          caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = Abs(normal*(pDef.XYZ()-p1));
-          if (defl > maxdef) maxdef = defl;
-          if (defl > defface)
-          {
-            if (isfirst) break;
-            if (caninsert)
-            {
-              // Record new vertex
-             aNbPnt++;              
-              nbLocat++;
-              Location3d.Bind(nbLocat,pDef);
-              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
-              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
-              InternalV.Append(InsVertex);
-            }
-            badTriangles.Append(triDom.Key());
-          }
-        }
-      }
-
-      if (!m3) // Not a boundary
-      {
-        // Check if this link was already processed
-        if (v1 < v3) { nf = v1; nl = v3; } else { nf = v3; nl = v1; }
-        if (theCouples.Add(MeshShape_Couple(nf,nl)))
-        {
-          // Check deflection on edge 2
-          mi2d = (xy3+xy1)*0.5;
-          caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = Abs(normal*(pDef.XYZ()-p2));
-          if (defl > maxdef) maxdef = defl;
-          if (defl > defface)
-          {
-            if (isfirst) break;
-            if (caninsert)
-            {
-              // Record new vertex
-              aNbPnt++;
-              nbLocat++;
-              Location3d.Bind(nbLocat,pDef);
-              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);
-              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
-              InternalV.Append(InsVertex);
-            }
-            badTriangles.Append(triDom.Key());
-          }
-        }
-      }
-
-      if (!m1) // Not a boundary
-      {
-        // Check if this link was already processed
-        if (v1 < v2) { nf = v1; nl = v2; } else { nf = v2; nl = v1; }
-        if (theCouples.Add(MeshShape_Couple(nf,nl)))
-        {
-          // Check deflection on edge 3
-          mi2d = (xy1+xy2)*0.5;
-          caro->D0(mi2d.X(), mi2d.Y(), pDef);
-          defl = Abs(normal*(pDef.XYZ()-p3));
-          if (defl > maxdef) maxdef = defl;
-          if (defl > defface)
-          {
-            if (isfirst) break;
-            if (caninsert)
-            {
-              // Record new vertex
-              aNbPnt++;
-              nbLocat++;
-              Location3d.Bind(nbLocat,pDef);
-              mi2d.SetCoord((mi2d.X()-myumin)/deltaX,(mi2d.Y()-myvmin)/deltaY);              
-              InsVertex.Initialize(mi2d,nbDomains,nbLocat,MeshDS_Free);
-              InternalV.Append(InsVertex);
-            }
-            badTriangles.Append(triDom.Key());
-          }
-        }
-      }
-    }
-
-    if (!isfirst && InternalV.Extent() > 0) 
-    {
-      BRepMesh_Array1OfVertexOfDelaun verttab(1, InternalV.Extent());
-      BRepMesh_ListIteratorOfListOfVertex itVer(InternalV);
-      Standard_Integer ipn = 1;
-      for (; itVer.More(); itVer.Next())
-        verttab(ipn++) = itVer.Value();
-      trigu.AddVertices(verttab);
-      nbInserted++;
-    }
-  }
-
-  return maxdef;
-}
-
-
-//=======================================================================
-//function : Plan
-//purpose  : 
-//=======================================================================
-static Standard_Boolean Plan(const MeshShape_SurfacePoint& P1,
-			     const MeshShape_SurfacePoint& P2,
-			     const MeshShape_SurfacePoint& P3,
-			     gp_XYZ&        eqPlan,
-			     Standard_Real& diPolr)
-{
-  gp_XYZ v1(P2.Coord()-P1.Coord());
-  if (v1.Modulus()>Precision::Confusion()) {
-    gp_XYZ v2(P3.Coord()-P2.Coord());
-    if (v2.Modulus()>Precision::Confusion()) {
-      gp_XYZ v3(P1.Coord()-P3.Coord());
-      if (v3.Modulus()>Precision::Confusion()) {
-	eqPlan=v1^v2;
-	if (eqPlan.Modulus()>0.) {
-	  eqPlan.Normalize();
-	  diPolr=P1.Coord()*eqPlan;
-	  return Standard_True;
-	}
-      }
-    }
-  }
-  diPolr=0.;
-  eqPlan.SetCoord(0., 0., 0.);
-  return Standard_False;
-}
-
-
-//=======================================================================
-//function : Uindex
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::Uindex(const MeshShape_SurfacePoint& P) const
-{
-//#ifndef DEB 
-  return (Standard_Integer) Round((P.UV().X()-myumin)/(myumax-myumin)*LIMITE_TRIANGULATION);
-//#else
-//  return Round((P.UV().X()-myumin)/(myumax-myumin)*LIMITE_TRIANGULATION);
-//#endif
-}
-
-
-//=======================================================================
-//function : Vindex
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::Vindex(const MeshShape_SurfacePoint& P) const
-{
-//#ifndef DEB
-  return (Standard_Integer) Round((P.UV().Y()-myvmin)/(myvmax-myvmin)*LIMITE_TRIANGULATION);
-//#else
-//  return Round((P.UV().Y()-myvmin)/(myvmax-myvmin)*LIMITE_TRIANGULATION);
-//#endif
-}
-
-
-//=======================================================================
-//function : Append
-//purpose  : 
-//=======================================================================
-
-void BRepMesh_FastDiscret::Append(MeshShape_ListOfSurfacePoint& pntsOnSurf,
-			      const MeshShape_SurfacePoint& P)
-{
-  Standard_Integer IU = Uindex(P), IV = Vindex(P);
-
-  // pas d append dans le cas d une frontiere.
-  if ((IU != 0)                    && 
-      (IV != 0)                    && 
-      (IU != LIMITE_TRIANGULATION) && 
-      (IV != LIMITE_TRIANGULATION)) {
-    MeshShape_Couple C(IU, IV);
-    if (!mymap.Contains(C)) {
-      mymap.Add(C);
-      pntsOnSurf.Append(P);
-    }
-  }
-}
-
-//=======================================================================
-//function : AddInShape
-//purpose  : 
-//=======================================================================
-void BRepMesh_FastDiscret::AddInShape(const TopoDS_Face&  face,
-				  const Standard_Real defface)
-{
-//  gp_Pnt Pt;
-
-  try{
-  //TColStd_MapIteratorOfMapOfInteger it;
-  MeshDS_MapOfInteger::Iterator it;
-
-  Standard_Integer e1, e2, e3, nTri;
-  Standard_Integer v1, v2, v3, iv1, iv2, iv3;
-  Standard_Integer i, index;
-  Standard_Boolean o1, o2, o3;
-  TopAbs_Orientation orFace = face.Orientation();
-
-  //const TColStd_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
-  const MeshDS_MapOfInteger& TriMap = structure->ElemOfDomain(nbDomains);
-  it.Initialize(TriMap);
-    
-  nTri = TriMap.Extent();
-
-  if (nTri != 0) {
-    
-    Poly_Array1OfTriangle Tri(1, nTri);
-    
-    i = 1;
-    
-    for (; it.More(); it.Next()) {
-      structure->GetElement(it.Key()).Edges(e1, e2, e3, o1, o2, o3);
-      
-      const BRepMesh_Edge& ve1=structure->GetLink(e1);
-      const BRepMesh_Edge& ve2=structure->GetLink(e2);
-      const BRepMesh_Edge& ve3=structure->GetLink(e3);
-      
-      if (o1) {
-	v1=ve1.FirstNode();
-      }
-      else {
-	v1=ve1.LastNode();
-      }
-      if (o2)
-      {
-        v2=ve2.FirstNode();
-	v3=ve2.LastNode();
-      }
-      else
-      {
-	v3=ve2.FirstNode();
-      	v2=ve2.LastNode();
-      }
-      
-      iv1 = myvemap.FindIndex(v1);
-      if (iv1 == 0) iv1 = myvemap.Add(v1);
-      iv2 = myvemap.FindIndex(v2);
-      if (iv2 == 0) iv2 = myvemap.Add(v2);
-      iv3 = myvemap.FindIndex(v3);
-      if (iv3 == 0) iv3 = myvemap.Add(v3);
-      
-      if (orFace == TopAbs_REVERSED) Tri(i++).Set(iv1, iv3, iv2);
-      else Tri(i++).Set(iv1, iv2, iv3);
-    }
-    
-    Standard_Integer nbVertices = myvemap.Extent();
-    Handle(Poly_Triangulation) T = new Poly_Triangulation(nbVertices, nTri, Standard_True);
-    Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
-    Trian = Tri;
-    TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
-    TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
-    
-    for (i = 1; i <= nbVertices; i++) {
-      index = myvemap.FindKey(i);
-      Nodes(i) = Pnt(index);
-      Nodes2d(i).SetXY(Vertex(index).Coord());
-    }
-    
-    T->Deflection(defface);
-    
-    // stockage de la triangulation dans la BRep.
-    BRep_Builder B;
-    TopLoc_Location l = face.Location();
-    if (!l.IsIdentity()) {
-      gp_Trsf tr = l.Transformation();
-      tr.Invert();
-      for (i = Nodes.Lower(); i <= Nodes.Upper(); i++) 
-	Nodes(i).Transform(tr);
-    }
-    B.UpdateFace(face, T);
-
-    // mise en place des polygones sur triangulation dans la face:
-    MeshShape_DataMapIteratorOfDataMapOfShapeListOfTransient It(internaledges);
-
-    for (; It.More(); It.Next()) {
-      if (It.Value().Extent() == 1) {
-	const Handle(Poly_PolygonOnTriangulation)& NOD = 
-	  *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().First()));
-	B.UpdateEdge(TopoDS::Edge(It.Key()), NOD, T, l);
-      }
-      else {
-	const Handle(Poly_PolygonOnTriangulation)& NOD1 = 
-	  *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().First()));
-	const Handle(Poly_PolygonOnTriangulation)& NOD2 = 
-	  *((Handle(Poly_PolygonOnTriangulation)*)&(It.Value().Last()));
-	B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, l);
-      }
-    }
-  }
- }
- catch(Standard_Failure)
- {
-   MESH_FAILURE(face);
- }
-}
-
-
-//=======================================================================
-//function : output
-//purpose  : 
-//=======================================================================
-Standard_Integer BRepMesh_FastDiscret::NbTriangles() const
-{
-  return structure->NbElements();
-}
-
-//=======================================================================
-//function : Triangle
-//purpose  : 
-//=======================================================================
-
-const BRepMesh_Triangle& BRepMesh_FastDiscret::Triangle
-  (const Standard_Integer Index) const
-{
-  return structure->GetElement(Index);
-}
-
-//=======================================================================
-//function : NbEdges
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::NbEdges() const
-{
-  return structure->NbLinks();
-}
-
-//=======================================================================
-//function : Edge
-//purpose  : 
-//=======================================================================
-
-const BRepMesh_Edge& BRepMesh_FastDiscret::Edge(const Standard_Integer Index) const
-{
-  return structure->GetLink(Index);
-}
-
-//=======================================================================
-//function : NbVertices
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::NbVertices() const
-{
-  return structure->NbNodes();
-}
-
-//=======================================================================
-//function : Vertex
-//purpose  : 
-//=======================================================================
-
-const BRepMesh_Vertex& BRepMesh_FastDiscret::Vertex
-  (const Standard_Integer Index) const
-{
-  return structure->GetNode(Index);
-}
-
-//=======================================================================
-//function : Pnt
-//purpose  : 
-//=======================================================================
-
-const gp_Pnt& BRepMesh_FastDiscret::Pnt(const Standard_Integer Index) const
-{
-  return Location3d(structure->GetNode(Index).Location3d());
-}
-
-//=======================================================================
-//function : NbPoint3d
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::NbPoint3d() const
-{
-  return nbLocat;
-}
-
-//=======================================================================
-//function : Point3d
-//purpose  : 
-//=======================================================================
-
-const gp_Pnt& BRepMesh_FastDiscret::Point3d(const Standard_Integer Index) const
-{
-  return Location3d(Index);
-}
-
-//=======================================================================
-//function : Normal
-//purpose  : 
-//=======================================================================
-
-void BRepMesh_FastDiscret::Normal(const Standard_Integer Index,
-				 gp_Pnt& Pon, gp_Dir& Nor) const
-{
-  const BRepMesh_Vertex& v=structure->GetNode(Index);
-  BRepMesh_GeomTool::Normal(new BRepAdaptor_HSurface(Domains(v.Domain())), 
-		   v.Coord().X(), v.Coord().Y(), Pon, Nor);
-}
-
-//=======================================================================
-//function : Result
-//purpose  : 
-//=======================================================================
-
-Handle(BRepMesh_DataStructureOfDelaun) BRepMesh_FastDiscret::Result() const
-{
-  return structure;
-}
-
-
-//=======================================================================
-//function : NbDomains
-//purpose  : 
-//=======================================================================
-
-Standard_Integer BRepMesh_FastDiscret::NbDomains() const 
-{ return nbDomains;}
-
-//=======================================================================
-//function : DomainFace
-//purpose  : 
-//=======================================================================
-
-const TopoDS_Face& BRepMesh_FastDiscret::DomainFace(const Standard_Integer Index) const 
-{ 
-  return Domains(Index); 
-}
-
-//=======================================================================
-//function : VerticesOfDomain
-//purpose  : 
-//=======================================================================
-
-void BRepMesh_FastDiscret::VerticesOfDomain(const Standard_Integer Index,
-					MeshDS_MapOfInteger&  Indices) const 
-{ 
-  Indices.Clear();
-  
-  // recup de la map des edges.
-  const MeshDS_MapOfInteger& edmap = structure->LinkOfDomain(Index);
-
-  // iterateur sur les edges.
-  //TColStd_MapIteratorOfMapOfInteger iter(edmap);
-  MeshDS_MapOfInteger::Iterator iter(edmap);
-  
-  Standard_Integer ind_edge;
-  for (iter.Reset(); iter.More(); iter.Next()) {
-    ind_edge = iter.Key();
-    const BRepMesh_Edge& Ed = Edge(ind_edge);
-    Indices.Add(Ed.FirstNode());
-    Indices.Add(Ed.LastNode());
-  }
-}
-
-
-
-//=======================================================================
-//function : EdgesOfDomain
-//purpose  : 
-//=======================================================================
-
-void BRepMesh_FastDiscret::EdgesOfDomain(const Standard_Integer Index,
-				     MeshDS_MapOfInteger&  Indices) const 
-{ 
-  Indices = structure->LinkOfDomain(Index);
-}
-
-
-//=======================================================================
-//function : TrianglesOfDomain
-//purpose  : 
-//=======================================================================
-
-void BRepMesh_FastDiscret::TrianglesOfDomain(const Standard_Integer Index,
-				     MeshDS_MapOfInteger&  Indices) const 
-{ 
-  Indices = structure->ElemOfDomain(Index);
-}
-
-
-
-//=======================================================================
-//function : GetDeflection
-//purpose  : 
-//=======================================================================
-
-Standard_Real BRepMesh_FastDiscret::GetDeflection() const
-{
-  return deflection;
-}
-
-
-//=======================================================================
-//function : GetAngle
-//purpose  : 
-//=======================================================================
-
-Standard_Real BRepMesh_FastDiscret::GetAngle() const
-{
-  return angle;
-}
-
-
-//=======================================================================
-//function : FindUV
-//purpose  : 
-//=======================================================================
-
-gp_XY BRepMesh_FastDiscret::FindUV(const TopoDS_Vertex&   V,
-			       const gp_Pnt2d&        XY, 
-			       const Standard_Integer ip,
-			       const Handle(BRepAdaptor_HSurface)& S,
-			       const Standard_Real mindist) 
-{
-  gp_XY theUV;
-  if (mylocation2d.IsBound(ip))
-  {
-    BRepMesh_ListOfXY& L = mylocation2d.ChangeFind(ip);
-    theUV = L.First();
-    if (L.Extent() != 1)
-    {
-      BRepMesh_ListIteratorOfListOfXY it(L);
-      it.Next();
-      Standard_Real dd, dmin = XY.Distance(gp_Pnt2d(theUV));
-      for (; it.More(); it.Next())
-      {
-	dd = XY.Distance(gp_Pnt2d(it.Value()));
-	if (dd < dmin)
-        {
-	  theUV = it.Value();
-	  dmin = dd;
-	}
-      }
-    }
-
-    const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(V), mindist);
-
-    const Standard_Real Utol2d = .5 * (S->LastUParameter() - S->FirstUParameter());
-    const Standard_Real Vtol2d = .5 * (S->LastVParameter() - S->FirstVParameter());
-
-    const gp_Pnt p1 = S->Value(theUV.X(), theUV.Y());
-    const gp_Pnt p2 = S->Value(XY.X(), XY.Y());
-
-    if (Abs(theUV.X() - XY.X()) > Utol2d ||
-        Abs(theUV.Y() - XY.Y()) > Vtol2d ||
-	!p1.IsEqual(p2, tol))
-    {
-      theUV = XY.Coord();
-      L.Append(theUV);
-    }
-  }
-  else
+  else
   {
     theUV = XY.Coord();
     BRepMesh_ListOfXY L;
@@ -3007,56 +1589,29 @@
   return myfacestate;
 }
 
-static Standard_Boolean GetVertexParameters(const TopoDS_Vertex& theVert, 
-                                         const TopoDS_Face& theFace,
-                                         gp_Pnt2d& thePoint)
+//=======================================================================
+//function : GetFaceAttribute
+//purpose  : 
+//=======================================================================
+
+Standard_Boolean BRepMesh_FastDiscret::GetFaceAttribute(const TopoDS_Face& theFace,
+                                                        Handle(BRepMesh_FaceAttribute)& fattrib) const
 {
-  TopLoc_Location L;
-  const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
-  L = L.Predivided(theVert.Location());
-  BRep_ListIteratorOfListOfPointRepresentation itpr
-    ((*((Handle(BRep_TVertex)*) &theVert.TShape()))->Points());
-  // On regarde dabord si il y des PointRepresentation (cas non Manifold)
-
-  while (itpr.More()) {
-    if (itpr.Value()->IsPointOnSurface(S,L)) {
-      thePoint.SetCoord(itpr.Value()->Parameter(),
-		     itpr.Value()->Parameter2());
-      return Standard_True;
-    }
-    itpr.Next();
+  if(mapattrib.IsBound(theFace))
+  {
+    fattrib = mapattrib(theFace);
+    return Standard_True;
   }
   return Standard_False;
 }
+
 //=======================================================================
-//function : Add
-//purpose  : method intended to addition internav vertices in triangulation.
+//function : RemoveFaceAttribute
+//purpose  : 
 //=======================================================================
 
-void BRepMesh_FastDiscret::Add(const TopoDS_Vertex&                theVert, 
-			       const TopoDS_Face&                  theFace, 
-			       const Handle(BRepAdaptor_HSurface)& thegFace)
-			       
+void BRepMesh_FastDiscret::RemoveFaceAttribute(const TopoDS_Face& theFace)
 {
-  const TopAbs_Orientation anOrient = theVert.Orientation();
-  gp_Pnt2d uvXY;
-  if( anOrient != TopAbs_INTERNAL || !GetVertexParameters(theVert,theFace,uvXY))
-    return;
-  Standard_Integer indVert =0;
-  if (vertices.IsBound(theVert))
-    indVert = vertices.Find(theVert);
-  else
-  {
-    nbLocat++;
-    Location3d.Bind(nbLocat, BRep_Tool::Pnt(theVert));
-    indVert = nbLocat;
-    vertices.Bind(theVert, indVert);
-  }
-  Standard_Real mindist = BRep_Tool::Tolerance(theVert);
-  // gp_Pnt2d uvXY = BRep_Tool::Parameters(theVert,theFace);
-  gp_XY anUV = FindUV(theVert, uvXY, indVert, thegFace, mindist);
-  BRepMesh_Vertex vf(anUV, nbDomains, indVert, MeshDS_Fixed);
-  Standard_Integer ivff = structure->AddNode(vf);
-  Standard_Integer isvf = myvemap.FindIndex(ivff);
-  if (isvf == 0) isvf = myvemap.Add(ivff);  
+  if(mapattrib.IsBound(theFace))
+    mapattrib.UnBind(theFace);
 }
