src/BRepMesh/BRepMesh_FastDiscretFace.cxx

   1 // File:        BRepMesh_FastDiscretFace.cxx
   2 // Created:
   3 // Author:      Ekaterina SMIRNOVA
   4 // Copyright: Open CASCADE SAS 2008
   5
   6 #include <BRepMesh_FastDiscretFace.ixx>
   7
   8 #include <Adaptor3d_IsoCurve.hxx>
   9 #include <BRepAdaptor_Curve.hxx>
  10 #include <BRepMesh_GeomTool.hxx>
  11 #include <BRepMesh_ListOfXY.hxx>
  12 #include <BRepMesh_Array1OfVertexOfDelaun.hxx>
  13 #include <BRepMesh_ListIteratorOfListOfVertex.hxx>
  14 #include <BRepMesh_ListIteratorOfListOfXY.hxx>
  15 #include <BRepMesh_PairOfPolygon.hxx>
  16 #include <BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon.hxx>
  17 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
  18 #include <BRepMesh_ClassifierPtr.hxx>
  19 #include <BRep_Builder.hxx>
  20 #include <BRep_PointRepresentation.hxx>
  21 #include <BRep_TVertex.hxx>
  22 #include <BRep_Tool.hxx>
  23 #include <Geom_Surface.hxx>
  24 #include <ElSLib.hxx>
  25 #include <Standard_ErrorHandler.hxx>
  26 #include <Standard_Failure.hxx>
  27 #include <Poly_PolygonOnTriangulation.hxx>
  28 #include <Poly_Triangulation.hxx>
  29 #include <TColStd_Array1OfReal.hxx>
  30 #include <TColStd_SequenceOfReal.hxx>
  31 #include <TColStd_Array1OfInteger.hxx>
  32 #include <TColStd_HArray1OfReal.hxx>
  33 #include <TopExp_Explorer.hxx>
  34 #include <TopoDS.hxx>
  35 #include <TopExp.hxx>
  36 #include <TColgp_Array1OfPnt2d.hxx>
  37 #include <NCollection_Map.hxx>
  38 #include <Geom_BSplineSurface.hxx>
  39 #include <GeomLib.hxx>
  40 #include <Bnd_Box2d.hxx>
  41
  42 static Standard_Real FUN_CalcAverageDUV(TColStd_Array1OfReal& P, const Standard_Integer PLen)
  43 {
  44   Standard_Integer i, j, n = 0;
  45   Standard_Real p, result = 0.;
  46
  47   for(i = 1; i <= PLen; i++)
  48   {
  49     // Sort
  50     for(j = i + 1; j <= PLen; j++)
  51     {
  52       if(P(i) > P(j))
  53       {
  54         p = P(i);
  55         P(i) = P(j);
  56         P(j) = p;
  57       }
  58     }
  59     // Accumulate
  60     if (i != 1)
  61     {
  62       p = Abs(P(i) - P(i-1));
  63       if(p > 1.e-7)
  64       {
  65         result += p;
  66         n++;
  67       }
  68     }
  69   }
  70   return (n? (result / (Standard_Real) n) : -1.);
  71 }
  72
  73 //=======================================================================
  74 //function : BRepMesh_FastDiscretFace
  75 //purpose  :
  76 //=======================================================================
  77 BRepMesh_FastDiscretFace::BRepMesh_FastDiscretFace
  78                           (const Standard_Real     theAngle,
  79                            const Standard_Boolean  theWithShare) :
  80   myAngle(theAngle), myWithShare(theWithShare), myNbLocat(0),
  81   myInternalVerticesMode(Standard_True)
  82 {
  83   myAllocator = new NCollection_IncAllocator(64000);
  84 }
  85
  86 //=======================================================================
  87 //function : Add(face)
  88 //purpose  :
  89 //=======================================================================
  90
  91 void BRepMesh_FastDiscretFace::Add(const TopoDS_Face&                    theFace,
  92                                    const Handle(BRepMesh_FaceAttribute)& theAttrib,
  93                                    const TopTools_DataMapOfShapeReal&    theMapDefle)
  94 {
  95 #ifndef DEB_MESH
  96   try
  97   {
  98     OCC_CATCH_SIGNALS
  99 #endif
 100     TopoDS_Face face = theFace;
 101     TopLoc_Location loc;
 102
 103     const Handle(Poly_Triangulation)& aFaceTrigu = BRep_Tool::Triangulation(face, loc);
 104     if ( aFaceTrigu.IsNull() )
 105       return;
 106
 107     myAttrib = theAttrib;
 108     face.Orientation(TopAbs_FORWARD);
 109     myStructure.Nullify();
 110     Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
 111     anAlloc->Reset(Standard_False);
 112     myStructure=new BRepMesh_DataStructureOfDelaun(anAlloc);
 113     BRepAdaptor_Surface  BS(face, Standard_False);
 114     Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
 115
 116     GeomAbs_SurfaceType thetype;
 117     thetype = BS.GetType();
 118
 119     TopAbs_Orientation orFace = face.Orientation();
 120
 121     if (!myWithShare)
 122       myVertices.Clear();
 123
 124     myListver.Clear();
 125     myVemap.Clear();
 126     myLocation2d.Clear();
 127     myInternaledges.Clear();
 128
 129     Standard_Integer i = 1;
 130     Standard_Integer ipn = 0;
 131
 132     TopoDS_Iterator exW(face);
 133     for (; exW.More(); exW.Next()) {
 134       const TopoDS_Shape& aWire = exW.Value();
 135       if (aWire.ShapeType() != TopAbs_WIRE)
 136         continue;
 137       TopoDS_Iterator ex(aWire);
 138       for(; ex.More(); ex.Next()) {
 139         const TopoDS_Edge& edge = TopoDS::Edge(ex.Value());
 140         if(edge.IsNull())
 141           continue;
 142
 143         RestoreStructureFromTriangulation(edge, face, gFace, aFaceTrigu, theMapDefle(edge), loc);
 144       }
 145     }
 146
 147     ////////////////////////////////////////////////////////////
 148     //add internal vertices after self-intersection check
 149     Standard_Integer nbVertices = 0;
 150     if(myInternalVerticesMode) {
 151       for(TopExp_Explorer ex(face,TopAbs_VERTEX ,TopAbs_EDGE); ex.More(); ex.Next()) {
 152         const TopoDS_Vertex& aVert = TopoDS::Vertex(ex.Current());
 153         Add(aVert,face,gFace);
 154       }
 155       nbVertices = myVemap.Extent();
 156     }
 157
 158     // essai de determination de la longueur vraie:
 159     // akm (bug OCC16) : We must calculate these measures in non-singular
 160     //     parts of face. Let`s try to compute average value of three
 161     //     (umin, (umin+umax)/2, umax), and respectively for v.
 162     //                 vvvvv
 163     //Standard_Real longu = 0.0, longv = 0.0; //, last , first;
 164     //gp_Pnt P11, P12, P21, P22, P31, P32;
 165
 166     Standard_Real umax = myAttrib->GetUMax();
 167     Standard_Real umin = myAttrib->GetUMin();
 168     Standard_Real vmax = myAttrib->GetVMax();
 169     Standard_Real vmin = myAttrib->GetVMin();
 170
 171     TColStd_Array1OfInteger tabvert_corr(1, nbVertices);
 172     gp_Pnt2d p2d;
 173
 174     // Check the necessity to fill the map of parameters
 175     const Standard_Boolean useUVParam = (thetype == GeomAbs_Torus ||
 176                                          thetype == GeomAbs_BezierSurface ||
 177                                          thetype == GeomAbs_BSplineSurface);
 178     myUParam.Clear();
 179     myVParam.Clear();
 180
 181     BRepMesh_IDMapOfNodeOfDataStructureOfDelaun aMoveNodes(myVemap.Extent());
 182
 183     for (i = 1; i <= myStructure->NbNodes(); i++)
 184     {
 185       const BRepMesh_Vertex& v = myStructure->GetNode(i);
 186       p2d = v.Coord();
 187       if (useUVParam) {
 188         myUParam.Add(p2d.X());
 189         myVParam.Add(p2d.Y());
 190       }
 191       gp_XY res;
 192       res.SetCoord((p2d.X()-(myAttrib->GetMinX()))/(myAttrib->GetDeltaX()),
 193                    (p2d.Y()-(myAttrib->GetMinY()))/(myAttrib->GetDeltaY()));
 194       BRepMesh_Vertex v_new(res,v.Location3d(),v.Movability());
 195       const BRepMesh_ListOfInteger& alist = myStructure->GetNodeList(i);
 196       aMoveNodes.Add(v_new,alist);
 197       tabvert_corr(i) = i;
 198     }
 199     myStructure->ReplaceNodes(aMoveNodes);
 200
 201     Standard_Boolean rajout;
 202
 203     BRepMesh_ClassifierPtr& classifier = theAttrib->GetClassifier();
 204
 205     switch (thetype)
 206     {
 207     case GeomAbs_Plane:
 208       rajout = !classifier->NaturalRestriction();
 209       break;
 210     case GeomAbs_Sphere:
 211     case GeomAbs_Torus:
 212       rajout = Standard_True;
 213       break;
 214     default:
 215       rajout = Standard_False;
 216     }
 217
 218     BRepMesh_Delaun trigu(myStructure, tabvert_corr, orFace==TopAbs_FORWARD);
 219
 220     //removed all free edges from triangulation
 221     Standard_Integer nbLinks = myStructure->NbNodes();
 222     for(i = 1; i <= nbLinks; i++)
 223     {
 224       if(myStructure->ElemConnectedTo(i).Extent() < 1)
 225       {
 226         BRepMesh_Edge& anEdge = (BRepMesh_Edge&)trigu.GetEdge(i);
 227         if(anEdge.Movability()==BRepMesh_Deleted)
 228           continue;
 229         anEdge.SetMovability(BRepMesh_Free);
 230         myStructure->RemoveLink(i);
 231       }
 232     }
 233
 234     Standard_Boolean isaline;
 235     isaline = ((umax-umin)<1.e-05) || ((vmax-vmin)<1.e-05);
 236
 237     Standard_Real aDef = -1;
 238     if (!isaline && myStructure->ElemOfDomain().Extent() > 0) {
 239       TColStd_ListOfInteger badTri, nulTri;
 240
 241       if(!rajout)
 242       {
 243         aDef = Control(gFace, theAttrib->GetDefFace(), myListver, badTri, nulTri, trigu, Standard_True);
 244         if( aDef > theAttrib->GetDefFace() || aDef < 0.)
 245           rajout = Standard_True;
 246       }
 247
 248       if(!rajout) {
 249         if(useUVParam) {
 250           if(BS.IsUClosed()) {
 251             if(myVParam.Extent() > 2) {
 252               rajout = Standard_True;
 253             }
 254           }
 255           if(BS.IsVClosed()) {
 256             if(myUParam.Extent() > 2) {
 257               rajout = Standard_True;
 258             }
 259           }
 260         }
 261       }
 262
 263       if(rajout){
 264         InternalVertices(gFace, myListver, theAttrib->GetDefFace(), classifier);
 265
 266         if (myListver.Extent() > 0) {
 267           BRepMesh_Array1OfVertexOfDelaun verttab(1, myListver.Extent());
 268           BRepMesh_ListIteratorOfListOfVertex itVer(myListver);
 269           ipn = 1;
 270           for (; itVer.More(); itVer.Next())
 271             verttab(ipn++) = itVer.Value();
 272           trigu.AddVertices(verttab);
 273         }
 274         //control internal points
 275         BRepMesh_ListOfVertex vvlist;
 276         aDef = Control(gFace, theAttrib->GetDefFace(), vvlist, badTri, nulTri, trigu, Standard_False);
 277         myListver.Append(vvlist);
 278       }
 279     }
 280
 281     //modify myStructure back
 282     aMoveNodes.Clear();
 283     Standard_Real deltaX = myAttrib->GetDeltaX();
 284     Standard_Real deltaY = myAttrib->GetDeltaY();
 285     for (i = 1; i <= myStructure->NbNodes(); i++)
 286     {
 287       const BRepMesh_Vertex& v = myStructure->GetNode(i);
 288       p2d = v.Coord();
 289       gp_XY res;
 290       res.SetCoord(p2d.X()*deltaX+umin,p2d.Y()*deltaY+vmin);
 291       BRepMesh_Vertex v_new(res,v.Location3d(),v.Movability());
 292       const BRepMesh_ListOfInteger& alist = myStructure->GetNodeList(i);
 293       aMoveNodes.Add(v_new,alist);
 294     }
 295     myStructure->ReplaceNodes(aMoveNodes);
 296
 297     AddInShape(face, (aDef < 0.0)? theAttrib->GetDefFace() : aDef);
 298 #ifndef DEB_MESH
 299   }
 300   catch(Standard_Failure)
 301   {
 302     BRep_Builder B;
 303     Handle(Poly_Triangulation) TNull;
 304     B.UpdateFace(theFace,TNull);
 305   }
 306 #endif // DEB_MESH
 307   myStructure.Nullify();
 308   myAttrib.Nullify();
 309 }
 310
 311 //=======================================================================
 312 //function : RestoreStructureFromTriangulation(edge)
 313 //purpose  : Restore structure of Delaun from triangulation on face
 314 //=======================================================================
 315 Standard_Boolean BRepMesh_FastDiscretFace::RestoreStructureFromTriangulation
 316                                (const TopoDS_Edge&                  theEdge,
 317                                 const TopoDS_Face&                  theFace,
 318                                 const Handle(BRepAdaptor_HSurface)& theSurf,
 319                                 const Handle(Poly_Triangulation)&   theTrigu,
 320                                 const Standard_Real                 theDefEdge,
 321                                 const TopLoc_Location&              theLoc)
 322 {
 323   // oan: changes for right restoring of triangulation data from face & edges
 324   Handle(Poly_PolygonOnTriangulation) Poly;
 325   Poly = BRep_Tool::PolygonOnTriangulation(theEdge, theTrigu, theLoc);
 326
 327   if (Poly.IsNull() || !Poly->HasParameters())
 328   {
 329     return Standard_False;
 330   }
 331
 332   // 2d vertex indices
 333   TopAbs_Orientation orEdge = theEdge.Orientation();
 334   // Get end points on 2d curve
 335   gp_Pnt2d uvFirst, uvLast;
 336   BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
 337
 338   // Get vertices
 339   TopoDS_Vertex pBegin, pEnd;
 340   TopExp::Vertices(theEdge,pBegin,pEnd);
 341
 342   const Standard_Boolean sameUV =
 343     uvFirst.IsEqual(uvLast, Precision::PConfusion());
 344
 345   const TColgp_Array1OfPnt2d&    UVNodes = theTrigu->UVNodes();
 346   const TColgp_Array1OfPnt&      Nodes   = theTrigu->Nodes();
 347   const TColStd_Array1OfInteger& Indices = Poly->Nodes();
 348
 349   const Standard_Integer nbnodes = Indices.Length();
 350   TColStd_Array1OfInteger NewNodes(1, nbnodes);
 351
 352   gp_Pnt  P3d;
 353   gp_XY   anUV;
 354
 355   // Process first vertex
 356   Standard_Integer ipf;
 357   if (myVertices.IsBound(pBegin))
 358   {
 359     ipf = myVertices.Find(pBegin);
 360   }
 361   else
 362   {
 363     if (sameUV && myVertices.IsBound(pEnd))
 364     {
 365       ipf = myVertices.Find(pEnd);
 366     }
 367     else
 368     {
 369       P3d = Nodes(Indices(1));
 370       if (!theLoc.IsIdentity())
 371         P3d.Transform(theLoc.Transformation());
 372       myNbLocat++;
 373       myLocation3d.Bind(myNbLocat, P3d);
 374       ipf = myNbLocat;
 375     }
 376     myVertices.Bind(pBegin,ipf);
 377   }
 378
 379  //Controle vertice tolerances
 380   gp_Pnt pFirst = theSurf->Value(uvFirst.X(), uvFirst.Y());
 381   gp_Pnt pLast  = theSurf->Value(uvLast.X(), uvLast.Y());
 382
 383   Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)),
 384                                     pLast.Distance(BRep_Tool::Pnt(pEnd)));
 385
 386   if (mindist < BRep_Tool::Tolerance(pBegin) ||
 387       mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
 388
 389   anUV = FindUV(pBegin, uvFirst, ipf, theSurf, mindist, myLocation2d);
 390   Standard_Integer iv1, isv1;
 391   BRepMesh_Vertex vf(anUV, ipf, BRepMesh_Frontier);
 392   iv1 = myStructure->AddNode(vf);
 393   isv1 = myVemap.FindIndex(iv1);
 394   if (isv1 == 0)
 395     isv1 = myVemap.Add(iv1);
 396   NewNodes(1) = isv1;
 397
 398   // Process last vertex
 399   Standard_Integer ipl, ivl;
 400   if (pEnd.IsSame(pBegin))
 401   {
 402     ipl = ipf;
 403   }
 404   else
 405   {
 406     if (myVertices.IsBound(pEnd))
 407     {
 408       ipl = myVertices.Find(pEnd);
 409     }
 410     else
 411     {
 412       if (sameUV)
 413       {
 414         ipl = ipf;
 415         ivl = iv1;
 416       }
 417       else
 418       {
 419         P3d = Nodes(Indices(nbnodes));
 420         if (!theLoc.IsIdentity())
 421           P3d.Transform(theLoc.Transformation());
 422         myNbLocat++;
 423         myLocation3d.Bind(myNbLocat, P3d);
 424         ipl = myNbLocat;
 425       }
 426       myVertices.Bind(pEnd,ipl);
 427     }
 428   }
 429
 430   anUV = FindUV(pEnd, uvLast, ipl, theSurf, mindist, myLocation2d);
 431   BRepMesh_Vertex vl(anUV, ipl, BRepMesh_Frontier);
 432
 433   Standard_Integer isvl;
 434   ivl = myStructure->AddNode(vl);
 435   isvl = myVemap.FindIndex(ivl);
 436   if (isvl == 0)
 437     isvl = myVemap.Add(ivl);
 438   NewNodes(nbnodes) = isvl;
 439
 440   BRepMesh_Vertex v;
 441
 442   Standard_Integer i;
 443   for (i = 2; i < nbnodes; i++)
 444   {
 445     // Record 3d point
 446     P3d = Nodes(Indices(i));
 447     if (!theLoc.IsIdentity())
 448       P3d.Transform(theLoc.Transformation());
 449     myNbLocat++;
 450     myLocation3d.Bind(myNbLocat, P3d);
 451
 452     // Record 2d point
 453     anUV = UVNodes(Indices(i)).Coord();
 454
 455     Standard_Integer iv2, isv;
 456     v.Initialize(anUV, myNbLocat, BRepMesh_Frontier);
 457     iv2 = myStructure->AddNode(v);
 458     isv = myVemap.FindIndex(iv2);
 459     if (isv == 0)
 460       isv = myVemap.Add(iv2);
 461     NewNodes(i) = isv;
 462
 463     //add links
 464     if (orEdge == TopAbs_FORWARD)
 465       myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
 466     else if (orEdge == TopAbs_REVERSED)
 467       myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
 468     else if (orEdge == TopAbs_INTERNAL)
 469       myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
 470     iv1 = iv2;
 471   }
 472
 473   // last point
 474   if (iv1 != ivl) {
 475     if (orEdge == TopAbs_FORWARD)
 476       myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
 477     else if (orEdge == TopAbs_REVERSED)
 478       myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
 479     else if (orEdge == TopAbs_INTERNAL)
 480       myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
 481   }
 482
 483   Handle(Poly_PolygonOnTriangulation) P1 =
 484         new Poly_PolygonOnTriangulation(NewNodes, Poly->Parameters()->Array1());
 485   P1->Deflection(theDefEdge);
 486   if (myInternaledges.IsBound(theEdge))
 487   {
 488     BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
 489     if (theEdge.Orientation() == TopAbs_REVERSED)
 490       pair.Append(P1);
 491     else
 492       pair.Prepend(P1);
 493   }
 494   else
 495   {
 496     BRepMesh_PairOfPolygon pair1;
 497     pair1.Append(P1);
 498     myInternaledges.Bind(theEdge, pair1);
 499   }
 500
 501   return Standard_True;
 502 }
 503
 504 //=======================================================================
 505 //function : InternalVertices
 506 //purpose  :
 507 //=======================================================================
 508 void BRepMesh_FastDiscretFace::InternalVertices(const Handle(BRepAdaptor_HSurface)& theCaro,
 509                                                 BRepMesh_ListOfVertex&              theInternalV,
 510                                                 const Standard_Real                 theDefFace,
 511                                                 const BRepMesh_ClassifierPtr&       theClassifier)
 512 {
 513   BRepMesh_Vertex newV;
 514   gp_Pnt2d p2d;
 515   gp_Pnt p3d;
 516
 517   // travail suivant le type de surface
 518   const BRepAdaptor_Surface& BS = *(BRepAdaptor_Surface*)&(theCaro->Surface());
 519   GeomAbs_SurfaceType thetype = theCaro->GetType();
 520
 521   Standard_Real umax = myAttrib->GetUMax();
 522   Standard_Real umin = myAttrib->GetUMin();
 523   Standard_Real vmax = myAttrib->GetVMax();
 524   Standard_Real vmin = myAttrib->GetVMin();
 525   Standard_Real deltaX = myAttrib->GetDeltaX();
 526   Standard_Real deltaY = myAttrib->GetDeltaY();
 527
 528   if (thetype == GeomAbs_Plane && !theClassifier->NaturalRestriction())
 529   {
 530     // rajout d`un seul point au milieu.
 531     const Standard_Real U = 0.5*(umin+umax);
 532     const Standard_Real V = 0.5*(vmin+vmax);
 533     if (theClassifier->Perform(gp_Pnt2d(U, V)) == TopAbs_IN)
 534     {
 535       // Record 3d point
 536       BRepMesh_GeomTool::D0(theCaro, U, V, p3d);
 537       myNbLocat++;
 538       myLocation3d.Bind(myNbLocat, p3d);
 539       // Record 2d point
 540       p2d.SetCoord((U-umin)/deltaX, (V-vmin)/deltaY);
 541       newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 542       theInternalV.Append(newV);
 543     }
 544   }
 545   else if (thetype == GeomAbs_Sphere)
 546   {
 547     gp_Sphere S = BS.Sphere();
 548     const Standard_Real R = S.Radius();
 549
 550     // Calculate parameters for iteration in V direction
 551     Standard_Real Dv = 1.0 - (theDefFace/R);
 552     if (Dv < 0.0) Dv = 0.0;
 553     Standard_Real oldDv = 2.0 * ACos (Dv);
 554     Dv  =  .7 * oldDv; //.7 ~= sqrt(2.) - Dv is hypotenuse of triangle when oldDv is legs
 555     const Standard_Real sv = vmax - vmin;
 556     Dv = sv/((Standard_Integer)(sv/Dv) + 1);
 557     const Standard_Real pasvmax = vmax-Dv*0.5;
 558
 559     //Du can be defined from relation: 2*r*Sin(Du/2) = 2*R*Sin(Dv/2), r = R*Cos(v)
 560     //here approximate relation r*Du = R*Dv is used
 561
 562     Standard_Real Du, pasu, pasv; //, ru;
 563     const Standard_Real su = umax-umin;
 564     Standard_Boolean Shift = Standard_False;
 565     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
 566     {
 567       // Calculate parameters for iteration in U direction
 568       // 1.-.365*pasv*pasv is simple approximation of Cos(pasv)
 569       // with condition that it gives ~.1 when pasv = pi/2
 570       Du = Dv/(1.-.365*pasv*pasv);
 571       Du = su/((Standard_Integer)(su/Du) + 1);
 572       Shift = !Shift;
 573       const Standard_Real d = (Shift)? Du*.5 : 0.;
 574       const Standard_Real pasumax = umax-Du*0.5 + d;
 575       for (pasu = umin + Du - d; pasu < pasumax; pasu += Du)
 576       {
 577         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 578         {
 579           // Record 3d point
 580 #ifdef DEB_MESH_CHRONO
 581           D0Internal++;
 582 #endif
 583           ElSLib::D0(pasu, pasv, S, p3d);
 584           myNbLocat++;
 585           myLocation3d.Bind(myNbLocat, p3d);
 586           // Record 2d point
 587           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 588           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 589           theInternalV.Append(newV);
 590         }
 591       }
 592     }
 593   }
 594   else if (thetype == GeomAbs_Cylinder)
 595   {
 596     gp_Cylinder S = BS.Cylinder();
 597     const Standard_Real R = S.Radius();
 598
 599     // Calculate parameters for iteration in U direction
 600     Standard_Real Du = 1.0 - (theDefFace/R);
 601     if (Du < 0.0) Du = 0.0;
 602     Du = 2.0 * ACos (Du);
 603     if (Du > myAngle) Du = myAngle;
 604     const Standard_Real su = umax - umin;
 605     const Standard_Integer nbU = (Standard_Integer)(su/Du);
 606     Du = su/(nbU+1);
 607
 608     // Calculate parameters for iteration in V direction
 609     const Standard_Real sv = vmax - vmin;
 610     Standard_Integer nbV = (Standard_Integer)( nbU*sv/(su*R) );
 611     nbV = Min(nbV, 100*nbU);
 612     Standard_Real Dv = sv/(nbV+1);
 613
 614     Standard_Real pasu, pasv, pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
 615     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
 616       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
 617         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 618         {
 619           // Record 3d point
 620           ElSLib::D0(pasu, pasv, S, p3d);
 621           myNbLocat++;
 622           myLocation3d.Bind(myNbLocat, p3d);
 623           // Record 2d point
 624           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 625           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 626           theInternalV.Append(newV);
 627         }
 628       }
 629     }
 630   }
 631   else if (thetype == GeomAbs_Cone)
 632   {
 633     Standard_Real R, RefR, SAng;
 634     gp_Cone C = BS.Cone();
 635     RefR = C.RefRadius();
 636     SAng = C.SemiAngle();
 637     R = Max(Abs(RefR+vmin*Sin(SAng)), Abs(RefR+vmax*Sin(SAng)));
 638     Standard_Real Du, Dv, pasu, pasv;
 639     Du = Max(1.0e0 - (theDefFace/R),0.0e0);
 640     Du  = (2.0 * ACos (Du));
 641     Standard_Integer nbU = (Standard_Integer) ( (umax-umin)/Du );
 642     Standard_Integer nbV = (Standard_Integer) ( nbU*(vmax-vmin)/((umax-umin)*R) );
 643     Du = (umax-umin)/(nbU+1);
 644     Dv = (vmax-vmin)/(nbV+1);
 645
 646     Standard_Real pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
 647     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
 648       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
 649         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 650         {
 651           // Record 3d point
 652           ElSLib::D0(pasu, pasv, C, p3d);
 653           myNbLocat++;
 654           myLocation3d.Bind(myNbLocat, p3d);
 655           // Record 2d point
 656           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 657           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 658           theInternalV.Append(newV);
 659         }
 660       }
 661     }
 662   }
 663   else if (thetype == GeomAbs_Torus)
 664   {
 665     gp_Torus T = BS.Torus();
 666
 667     Standard_Boolean insert;
 668     Standard_Integer i, j, ParamULength, ParamVLength;
 669     Standard_Real pp, pasu, pasv;
 670     Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
 671
 672     TColStd_SequenceOfReal ParamU, ParamV;
 673
 674     Standard_Real Du, Dv;//, pasu, pasv;
 675     Dv = Max(1.0e0 - (theDefFace/r),0.0e0) ;
 676     Standard_Real oldDv = 2.0 * ACos (Dv);
 677     oldDv = Min(oldDv, myAngle);
 678     Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
 679     Dv = oldDv;
 680
 681     Standard_Integer nbV = Max((Standard_Integer)((vmax-vmin)/Dv), 2);
 682     Dv = (vmax-vmin)/(nbV+1);
 683     Standard_Real ru = R + r;
 684     if (ru > 1.e-16)
 685     {
 686       Du  = 2.0 * ACos(Max(1.0 - (theDefFace/ru),0.0));
 687       if (myAngle < Du) Du = myAngle;
 688       Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
 689       if(aa < gp::Resolution())
 690         return;
 691       Du *= Min(oldDv, Du) / aa;
 692     }
 693     else Du = Dv;
 694
 695     Standard_Integer nbU = Max((Standard_Integer)((umax-umin)/Du), 2);
 696     nbU = Max(nbU , (int)(nbV*(umax-umin)*R/((vmax-vmin)*r)/5.));
 697     Du = (umax-umin)/(nbU+1);
 698
 699     if (R < r)
 700     {
 701       // comme on recupere les points des edges.
 702       // dans ce cas, les points ne sont pas representatifs.
 703
 704       //-- On choisit DeltaX et DeltaY de facon a ce qu on ne saute pas
 705       //-- de points sur la grille
 706       for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
 707     }//R<r
 708     else //U if R > r
 709     {
 710       //--ofv: U
 711       // Number of mapped U parameters
 712       const Standard_Integer LenU = myUParam.Extent();
 713       // Fill array of U parameters
 714       TColStd_Array1OfReal Up(1,LenU);
 715       for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
 716
 717       // Calculate DU, sort array of parameters
 718       Standard_Real aDU = FUN_CalcAverageDUV(Up,LenU);
 719       aDU = Max(aDU, Abs(umax -  umin) / (Standard_Real) nbU / 2.);
 720       Standard_Real dUstd = Abs(umax -  umin) / (Standard_Real) LenU;
 721       if (aDU > dUstd) dUstd = aDU;
 722       // Add U parameters
 723       for (j = 1; j <= LenU; j++)
 724       {
 725         pp = Up(j);
 726         insert = Standard_True;
 727         ParamULength = ParamU.Length();
 728         for (i = 1; i <= ParamULength && insert; i++)
 729         {
 730           insert = (Abs(ParamU.Value(i)-pp) > (0.5*dUstd));
 731         }
 732         if (insert) ParamU.Append(pp);
 733       }
 734     }
 735
 736     //--ofv: V
 737     // Number of mapped V parameters
 738     const Standard_Integer LenV = myVParam.Extent();
 739     // Fill array of V parameters
 740     TColStd_Array1OfReal Vp(1,LenV);
 741     for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
 742     // Calculate DV, sort array of parameters
 743     Standard_Real aDV = FUN_CalcAverageDUV(Vp,LenV);
 744     aDV = Max(aDV, Abs(vmax -  vmin) / (Standard_Real) nbV / 2.);
 745
 746     Standard_Real dVstd = Abs(vmax -  vmin) / (Standard_Real) LenV;
 747     if (aDV > dVstd) dVstd = aDV;
 748     // Add V parameters
 749     for (j = 1; j <= LenV; j++)
 750     {
 751       pp = Vp(j);
 752
 753       insert = Standard_True;
 754       ParamVLength = ParamV.Length();
 755       for (i = 1; i <= ParamVLength && insert; i++)
 756       {
 757         insert = (Abs(ParamV.Value(i)-pp) > (dVstd*2./3.));
 758       }
 759       if (insert) ParamV.Append(pp);
 760     }
 761
 762     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
 763     Standard_Real uminnew = umin+deltaY*0.1;
 764     Standard_Real vminnew = vmin+deltaX*0.1;
 765     Standard_Real umaxnew = umax-deltaY*0.1;
 766     Standard_Real vmaxnew = vmax-deltaX*0.1;
 767
 768     for (i = 1; i <= Lu; i++)
 769     {
 770       pasu = ParamU.Value(i);
 771       if (pasu >= uminnew && pasu < umaxnew)
 772       {
 773         for (j = 1; j <= Lv; j++)
 774         {
 775           pasv = ParamV.Value(j);
 776           if (pasv >= vminnew && pasv < vmaxnew)
 777           {
 778             if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 779             {
 780               // Record 3d point
 781               ElSLib::D0(pasu, pasv, T, p3d);
 782               myNbLocat++;
 783               myLocation3d.Bind(myNbLocat, p3d);
 784               // Record 2d point
 785               p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 786               newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 787               theInternalV.Append(newV);
 788             }
 789           }
 790         }
 791       }
 792     }
 793   }
 794   else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
 795   {
 796     Standard_Integer i, j;
 797
 798     Handle(TColStd_HArray1OfReal) anUKnots, anVKnots;
 799     Standard_Integer aNbUNkots,aNbVNkots;
 800     Handle(Geom_Surface) B;
 801     if (thetype == GeomAbs_BezierSurface) {
 802       anUKnots = new TColStd_HArray1OfReal(1,2);
 803       anVKnots = new TColStd_HArray1OfReal(1,2);
 804       anUKnots->SetValue(1,0.);
 805       anUKnots->SetValue(2,1.);
 806       anVKnots->SetValue(1,0.);
 807       anVKnots->SetValue(2,1.);
 808       aNbUNkots = 2;
 809       aNbVNkots = 2;
 810       B = BS.Bezier();
 811     }
 812     else {
 813       Handle(Geom_BSplineSurface) aSurf = BS.BSpline();
 814       B = aSurf;
 815       aNbUNkots = aSurf->NbUKnots();
 816       aNbVNkots = aSurf->NbVKnots();
 817       anUKnots = new TColStd_HArray1OfReal(1,aNbUNkots);
 818       anVKnots = new TColStd_HArray1OfReal(1,aNbVNkots);
 819
 820       aSurf->UKnots(anUKnots->ChangeArray1());
 821       aSurf->VKnots(anVKnots->ChangeArray1());
 822     }
 823
 824     Standard_Real ddu = theCaro->UResolution(theDefFace)*5.;
 825     Standard_Real aDUmin = (umax-umin)/5.;
 826     if(ddu > aDUmin)
 827       ddu = aDUmin;
 828
 829     // Sort sequence of U parameters
 830     TColStd_SequenceOfReal ParamU;
 831     Standard_Integer anUdegree = theCaro->UDegree();
 832
 833     Standard_Real aU1 = anUKnots->Value(1);
 834     for( i = 1; i < aNbUNkots; i++)
 835     {
 836       Standard_Real aStart = anUKnots->Value(i);
 837       Standard_Real aEnd = anUKnots->Value(i+1);
 838       Standard_Integer aNbPoints = anUdegree-1;
 839       Standard_Real aStep = (aEnd-aStart)/(aNbPoints+1);
 840
 841       for( Standard_Integer anInd = 0; anInd<= aNbPoints; anInd++)
 842       {
 843         Standard_Real aU2 = aStart+anInd*aStep;
 844         if((aU2 - aU1) > ddu)
 845         {
 846           ParamU.Append(aU2);
 847           aU1 = aU2;
 848         }
 849       }
 850     }
 851     ParamU.Append(anUKnots->Value(aNbUNkots));
 852     Standard_Integer ParamULength = ParamU.Length();
 853
 854     Standard_Real ddv = theCaro->VResolution(theDefFace)*5.;
 855     Standard_Real aDVmin = (vmax-vmin)/5.;
 856     if(ddv > aDVmin)
 857       ddv = aDVmin;
 858     // Sort sequence of V parameters
 859     TColStd_SequenceOfReal ParamV;
 860     Standard_Integer anVdegree = theCaro->VDegree();
 861
 862     Standard_Real aV1 = anVKnots->Value(1);
 863     for( i = 1; i < aNbVNkots; i++)
 864     {
 865       Standard_Real aStart = anVKnots->Value(i);
 866       Standard_Real aEnd = anVKnots->Value(i+1);
 867       Standard_Integer aNbPoints = anVdegree-1;
 868       Standard_Real aStep = (aEnd-aStart)/(aNbPoints+1);
 869
 870       for( Standard_Integer anInd = 0; anInd<= aNbPoints; anInd++)
 871       {
 872         Standard_Real aV2 = aStart+anInd*aStep;
 873         if ((aV2 - aV1) > ddv)
 874         {
 875           ParamV.Append(aV2);
 876           aV1 = aV2;
 877         }
 878       }
 879     }
 880     ParamV.Append(anVKnots->Value(aNbVNkots));
 881     Standard_Integer ParamVLength = ParamV.Length();
 882
 883     // controle des isos U et insertion eventuelle:
 884
 885     gp_Pnt P1, P2, PControl;
 886     Standard_Real u, v, dist, V1, V2, U1, U2;
 887
 888     // precision for compare square distances
 889     double dPreci = Precision::Confusion()*Precision::Confusion();
 890
 891     // Insert V parameters by deflection criterion
 892     for (i = 1; i <= ParamULength; i++) {
 893       Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
 894       V1 = ParamV.Value(1);
 895       P1 = IsoU->Value(V1);
 896       for (j = 2; j <= ParamVLength;) {
 897         V2 = ParamV.Value(j);
 898         P2 = IsoU->Value(V2);
 899         v = 0.5*(V1+V2);
 900         PControl = IsoU->Value(v);
 901         // 23.03.2010 skl for OCC21645 - change precision for comparison
 902         if( P1.SquareDistance(P2) > dPreci ) {
 903           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
 904           dist = L.Distance(PControl);
 905         }
 906         else {
 907           dist = P1.Distance(PControl);
 908         }
 909         if (dist > theDefFace) {
 910           // insertion
 911           ParamV.InsertBefore(j, v);
 912           ParamVLength++;
 913         }
 914         else {
 915           V1 = V2;
 916           P1 = P2;
 917           j++;
 918         }
 919       }
 920     }
 921
 922     for (i = 2; i < ParamVLength; i++) {
 923       v = ParamV.Value(i);
 924       Handle(Geom_Curve) IsoV = B->VIso(v);
 925       U1 = ParamU.Value(1);
 926       P1 = IsoV->Value(U1);
 927       for (j = 2; j <= ParamULength;) {
 928         U2 = ParamU.Value(j);
 929         P2 = IsoV->Value(U2);
 930         u = 0.5*(U1+U2);
 931         PControl = IsoV->Value(u);
 932         // 23.03.2010 skl for OCC21645 - change precision for comparison
 933         if( P1.SquareDistance(P2) > dPreci ) {
 934           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
 935           dist = L.Distance(PControl);
 936         }
 937         else {
 938           dist = P1.Distance(PControl);
 939         }
 940         if (dist > theDefFace) {
 941           // insertion
 942           ParamU.InsertBefore(j, u);
 943           ParamULength++;
 944         }
 945         else {
 946           //check normal
 947           U1 = U2;
 948           P1 = P2;
 949           if (j < ParamULength) {
 950             // Classify intersection point
 951             if (theClassifier->Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
 952             {
 953               // Record 3d point
 954               myNbLocat++;
 955               myLocation3d.Bind(myNbLocat, P1);
 956               // Record 2d point
 957               p2d.SetCoord((U1-umin)/deltaX, (v-vmin)/deltaY);
 958               newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 959               theInternalV.Append(newV);
 960             }
 961           }
 962           j++;
 963         }
 964       }
 965     }
 966   }
 967   else {
 968     const Standard_Real anAngle = 0.35;
 969
 970     Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
 971     Adaptor3d_IsoCurve tabu[10], tabv[10];
 972
 973     TColStd_SequenceOfReal ParamU, ParamV;
 974     Standard_Real u, v, du, dv;
 975     Standard_Integer iu, iv;
 976     Standard_Real f, l;
 977
 978     du = (umax-umin) / (nbpointsU+1); dv = (vmax-vmin) / (nbpointsV+1);
 979
 980     for (iu = 1; iu <= nbpointsU; iu++) {
 981       u = umin + iu*du;
 982       tabu[iu-1].Load(theCaro);
 983       tabu[iu-1].Load(GeomAbs_IsoU, u);
 984     }
 985
 986     for (iv = 1; iv <= nbpointsV; iv++) {
 987       v = vmin + iv*dv;
 988       tabv[iv-1].Load(theCaro);
 989       tabv[iv-1].Load(GeomAbs_IsoV, v);
 990     }
 991
 992     Standard_Integer imax = 1, MaxV = 0;
 993
 994     GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU];
 995
 996     for (i = 0; i <= nbpointsU-1; i++) {
 997       f = Max(vmin, tabu[i].FirstParameter());
 998       l = Min(vmax, tabu[i].LastParameter());
 999       GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, anAngle, 0.7*theDefFace, 2);
1000       tabGU[i] = theDeflection;
1001       if (tabGU[i].NbPoints() > MaxV) {
1002         MaxV = tabGU[i].NbPoints();
1003         imax = i;
1004       }
1005     }
1006
1007     // recuperation du tableau de parametres V:
1008     Standard_Integer NV = tabGU[imax].NbPoints();
1009     for (i = 1; i <= NV; i++) {
1010       ParamV.Append(tabGU[imax].Parameter(i));
1011     }
1012     delete [] tabGU;
1013
1014     imax = 1;
1015     Standard_Integer MaxU = 0;
1016
1017     GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV];
1018
1019     for (i = 0; i <= nbpointsV-1; i++) {
1020       f = Max(umin, tabv[i].FirstParameter());
1021       l = Min(umax, tabv[i].LastParameter());
1022       GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, anAngle, 0.7*theDefFace, 2);
1023       tabGV[i] = thedeflection2;
1024       if (tabGV[i].NbPoints() > MaxU) {
1025         MaxU = tabGV[i].NbPoints();
1026         imax = i;
1027       }
1028     }
1029
1030     // recuperation du tableau de parametres U:
1031     Standard_Integer NU = tabGV[imax].NbPoints();
1032     for (i = 1; i <= NU; i++) {
1033       ParamU.Append(tabGV[imax].Parameter(i));
1034     }
1035     delete [] tabGV;
1036
1037     if (ParamU.Length() == 2) {
1038       ParamU.InsertAfter(1, (umax+umin)*0.5);
1039     }
1040     if (ParamV.Length() == 2) {
1041       ParamV.InsertAfter(1, (vmax+vmin)*0.5);
1042     }
1043
1044     TColStd_SequenceOfReal InsertV, InsertU;
1045     gp_Pnt P1;
1046
1047     Adaptor3d_IsoCurve IsoV;
1048     IsoV.Load(theCaro);
1049
1050     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
1051
1052     for (i = 2; i < Lv; i++) {
1053       v = ParamV.Value(i);
1054       IsoV.Load(GeomAbs_IsoV, v);
1055       for (j = 2; j < Lu; j++) {
1056         u = ParamU.Value(j);
1057         if (theClassifier->Perform(gp_Pnt2d(u, v)) == TopAbs_IN)
1058         {
1059           // Record 3d point
1060           P1 = IsoV.Value(u);
1061           myNbLocat++;
1062           myLocation3d.Bind(myNbLocat, P1);
1063           // Record 2d point
1064           p2d.SetCoord((u-umin)/deltaX, (v-vmin)/deltaY);
1065           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
1066           theInternalV.Append(newV);
1067         }
1068       }
1069     }
1070   }
1071 #ifdef DEB_MESH_CHRONO
1072   chInternal.Stop();
1073 #endif
1074 }
1075
1076 /**
1077 *  Internal class Couple, moved from MeshData package
1078 */
1079
1080 class BRepMesh_Couple
1081 {
1082  public:
1083   BRepMesh_Couple() { myI1 = myI2 = 0; }
1084   BRepMesh_Couple(const Standard_Integer I1,
1085                   const Standard_Integer I2)
1086   { myI1 = I1; myI2 = I2; }
1087
1088   Standard_Integer myI1;
1089   Standard_Integer myI2;
1090 };
1091
1092 inline Standard_Boolean IsEqual(const BRepMesh_Couple& one,
1093                                 const BRepMesh_Couple& other)
1094 {
1095   if (one.myI1 == other.myI1 &&
1096       one.myI2 == other.myI2) return Standard_True;
1097   else return Standard_False;
1098 }
1099
1100 inline Standard_Integer HashCode(const BRepMesh_Couple& one,
1101                                  const Standard_Integer Upper)
1102 {
1103   return ::HashCode((one.myI1+one.myI2), Upper);
1104 }
1105
1106 typedef NCollection_Map<BRepMesh_Couple> BRepMesh_MapOfCouple;
1107
1108 //=======================================================================
1109 //function : Control
1110 //purpose  :
1111 //=======================================================================
1112 Standard_Real BRepMesh_FastDiscretFace::Control(const Handle(BRepAdaptor_HSurface)& theCaro,
1113                                                 const Standard_Real                 theDefFace,
1114                                                 BRepMesh_ListOfVertex&              theInternalV,
1115                                                 TColStd_ListOfInteger&              theBadTriangles,
1116                                                 TColStd_ListOfInteger&              theNulTriangles,
1117                                                 BRepMesh_Delaun&                    theTrigu,
1118                                                 const Standard_Boolean              theIsFirst)
1119 {
1120   //IMPORTANT: Constants used in calculations
1121   const Standard_Real MinimalArea2d = 1.e-9;
1122   const Standard_Real MinimalSqLength3d = 1.e-12;
1123   const Standard_Real aDef2 = theDefFace*theDefFace;
1124
1125   // Define the number of iterations
1126   Standard_Integer myNbIterations = 11;
1127   const Standard_Integer nbPasses = (theIsFirst? 1 : myNbIterations);
1128
1129   // Initialize stop condition
1130   Standard_Boolean allDegenerated = Standard_False;
1131   Standard_Integer nbInserted = 1;
1132
1133   // Create map of links to skip already processed
1134   Standard_Integer nbtriangles = myStructure->ElemOfDomain().Extent();
1135   if (nbtriangles <= 0) return -1.0;
1136   BRepMesh_MapOfCouple theCouples(3*nbtriangles);
1137
1138   gp_XY mi2d;
1139   gp_XYZ vecEd1, vecEd2, vecEd3;
1140   gp_Pnt pDef;
1141   Standard_Real dv = 0, defl = 0, maxdef = -1;
1142   Standard_Integer pass = 1, nf = 0, nl = 0;
1143   BRepMesh_Vertex InsVertex;
1144   Standard_Boolean caninsert;
1145
1146   Standard_Real sqdefface = theDefFace * theDefFace;
1147   Standard_Real ddu = theCaro->UResolution(theDefFace);
1148   Standard_Real ddv = theCaro->VResolution(theDefFace);
1149
1150   GeomAbs_SurfaceType thetype = theCaro->GetType();
1151   Handle(Geom_Surface) BSpl;
1152   Standard_Boolean isSpline = Standard_False;
1153   if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
1154   {
1155     isSpline = Standard_True;
1156     if (thetype == GeomAbs_BezierSurface)
1157       BSpl = theCaro->Bezier();
1158     else
1159       BSpl = theCaro->BSpline();
1160   }
1161
1162   NCollection_DataMap<Standard_Integer,gp_Dir> aNorMap;
1163   NCollection_DataMap<Standard_Integer,Standard_Integer> aStatMap;
1164
1165   // Perform refinement passes
1166   for (; pass <= nbPasses && nbInserted && !allDegenerated; pass++)
1167   {
1168     theInternalV.Clear();
1169     theBadTriangles.Clear();
1170
1171     // Reset stop condition
1172     allDegenerated = Standard_True;
1173     nbInserted = 0;
1174     maxdef = -1.0;
1175
1176     // Do not insert nodes in last pass in non-SharedMode
1177     caninsert = (myWithShare || pass < nbPasses);
1178
1179     // Read mesh size
1180     nbtriangles = myStructure->ElemOfDomain().Extent();
1181     if (nbtriangles <= 0) break;
1182
1183     // Iterate on current triangles
1184     BRepMesh_MapOfInteger::Iterator triDom;
1185     const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1186     triDom.Initialize(TriMap);
1187     Standard_Integer aNbPnt = 0;
1188     Standard_Real umin = myAttrib->GetUMin();
1189     Standard_Real vmin = myAttrib->GetVMin();
1190     Standard_Real deltaX = myAttrib->GetDeltaX();
1191     Standard_Real deltaY = myAttrib->GetDeltaY();
1192     for (; triDom.More(); triDom.Next())
1193     {
1194       Standard_Integer TriId = triDom.Key();
1195       const BRepMesh_Triangle& curTri=Triangle(TriId);
1196       if (curTri.Movability()==BRepMesh_Deleted) continue;
1197
1198       Standard_Boolean o1, o2, o3;
1199       Standard_Integer v1 = 0, v2 = 0, v3 = 0, e1 = 0, e2 = 0, e3 = 0;
1200       curTri.Edges(e1, e2, e3, o1, o2, o3);
1201
1202       const BRepMesh_Edge& edg1=Edge(e1);
1203       const BRepMesh_Edge& edg2=Edge(e2);
1204       const BRepMesh_Edge& edg3=Edge(e3);
1205
1206       Standard_Boolean m1 = (edg1.Movability() == BRepMesh_Frontier);
1207       Standard_Boolean m2 = (edg2.Movability() == BRepMesh_Frontier);
1208       Standard_Boolean m3 = (edg3.Movability() == BRepMesh_Frontier);
1209       if (o1) {
1210         v1=edg1.FirstNode();
1211         v2=edg1.LastNode();
1212       }
1213       else {
1214         v1=edg1.LastNode();
1215         v2=edg1.FirstNode();
1216       }
1217       if (o2)
1218         v3=edg2.LastNode();
1219       else
1220         v3=edg2.FirstNode();
1221
1222       const BRepMesh_Vertex& vert1=Vertex(v1);
1223       const BRepMesh_Vertex& vert2=Vertex(v2);
1224       const BRepMesh_Vertex& vert3=Vertex(v3);
1225
1226       const gp_XYZ& p1=myLocation3d(vert1.Location3d()).Coord();
1227       const gp_XYZ& p2=myLocation3d(vert2.Location3d()).Coord();
1228       const gp_XYZ& p3=myLocation3d(vert3.Location3d()).Coord();
1229
1230       vecEd1 = p2 - p1;
1231       vecEd2 = p3 - p2;
1232       vecEd3 = p1 - p3;
1233
1234       // Check for degenerated triangle
1235       if (vecEd1.SquareModulus() < MinimalSqLength3d ||
1236           vecEd2.SquareModulus() < MinimalSqLength3d ||
1237           vecEd3.SquareModulus() < MinimalSqLength3d)
1238       {
1239         theNulTriangles.Append(TriId);
1240         continue;
1241       }
1242
1243       allDegenerated = Standard_False;
1244
1245       gp_XY xy1(vert1.Coord().X()*deltaX+umin,vert1.Coord().Y()*deltaY+vmin);
1246       gp_XY xy2(vert2.Coord().X()*deltaX+umin,vert2.Coord().Y()*deltaY+vmin);
1247       gp_XY xy3(vert3.Coord().X()*deltaX+umin,vert3.Coord().Y()*deltaY+vmin);
1248
1249       // Check triangle area in 2d
1250       if (Abs((xy2-xy1)^(xy3-xy1)) < MinimalArea2d)
1251       {
1252         theNulTriangles.Append(TriId);
1253         continue;
1254       }
1255
1256       // Check triangle normal
1257       gp_XYZ normal(vecEd1^vecEd2);
1258       dv = normal.Modulus();
1259       if (dv < Precision::Confusion())
1260       {
1261         theNulTriangles.Append(TriId);
1262         continue;
1263       }
1264       normal /= dv;
1265
1266       // Check deflection on triangle
1267       mi2d = (xy1+xy2+xy3)/3.0;
1268       theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1269       defl = pDef.SquareDistance((p1+p2+p3)/3.);
1270       if (defl > maxdef) maxdef = defl;
1271       if (defl > sqdefface)
1272       {
1273         if (theIsFirst) break;
1274         if (caninsert)
1275         {
1276           // Record new vertex
1277           aNbPnt++;
1278           myNbLocat++;
1279           myLocation3d.Bind(myNbLocat,pDef);
1280           mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1281           InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1282           theInternalV.Append(InsVertex);
1283         }
1284         theBadTriangles.Append(TriId);
1285       }
1286
1287       if (!m2) // Not a boundary
1288       {
1289         // Check if this link was already processed
1290         if (v2 < v3) { nf = v2; nl = v3; } else { nf = v3; nl = v2; }
1291         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1292         {
1293           // Check deflection on edge 1
1294           mi2d = (xy2+xy3)*0.5;
1295           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1296           defl = pDef.SquareDistance((p2+p3)/2.);
1297           if (defl > maxdef) maxdef = defl;
1298           if (defl > sqdefface)
1299           {
1300             if (theIsFirst) break;
1301             if (caninsert)
1302             {
1303               // Record new vertex
1304               aNbPnt++;
1305               myNbLocat++;
1306               myLocation3d.Bind(myNbLocat,pDef);
1307               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1308               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1309               theInternalV.Append(InsVertex);
1310             }
1311             theBadTriangles.Append(TriId);
1312           }
1313         }
1314       }
1315
1316       if (!m3) // Not a boundary
1317       {
1318         // Check if this link was already processed
1319         if (v1 < v3) { nf = v1; nl = v3; } else { nf = v3; nl = v1; }
1320         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1321         {
1322           // Check deflection on edge 2
1323           mi2d = (xy3+xy1)*0.5;
1324           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1325           defl = pDef.SquareDistance((p1+p3)/2.);
1326           if (defl > maxdef) maxdef = defl;
1327           if (defl > sqdefface)
1328           {
1329             if (theIsFirst) break;
1330             if (caninsert)
1331             {
1332               // Record new vertex
1333               aNbPnt++;
1334               myNbLocat++;
1335               myLocation3d.Bind(myNbLocat,pDef);
1336               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1337               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1338               theInternalV.Append(InsVertex);
1339             }
1340             theBadTriangles.Append(TriId);
1341           }
1342         }
1343       }
1344
1345       if (!m1) // Not a boundary
1346       {
1347         // Check if this link was already processed
1348         if (v1 < v2) { nf = v1; nl = v2; } else { nf = v2; nl = v1; }
1349         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1350         {
1351           // Check deflection on edge 3
1352           mi2d = (xy1+xy2)*0.5;
1353           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1354           defl = pDef.SquareDistance((p1+p2)/2.);
1355           if (defl > maxdef) maxdef = defl;
1356           if (defl > sqdefface)
1357           {
1358             if (theIsFirst) break;
1359             if (caninsert)
1360             {
1361               // Record new vertex
1362               aNbPnt++;
1363               myNbLocat++;
1364               myLocation3d.Bind(myNbLocat,pDef);
1365               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1366               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1367               theInternalV.Append(InsVertex);
1368             }
1369             theBadTriangles.Append(TriId);
1370           }
1371         }
1372       }
1373
1374       //check normal
1375       if(isSpline && !BSpl.IsNull() && (theBadTriangles.IsEmpty() || theBadTriangles.Last() != TriId))
1376       {
1377         gp_Dir N1(0,0,1), N2(0,0,1), N3(0,0,1);
1378         Standard_Integer aSt1, aSt2, aSt3;
1379         if(aNorMap.IsBound(v1)) {
1380           aSt1 = aStatMap.Find(v1);
1381           N1 =aNorMap.Find(v1);
1382         }
1383         else {
1384           aSt1 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy1), Precision::Confusion(), N1);
1385           aStatMap.Bind(v1,aSt1);
1386           aNorMap.Bind(v1,N1);
1387         }
1388
1389         if(aNorMap.IsBound(v2)) {
1390           aSt2 = aStatMap.Find(v2);
1391           N2 = aNorMap.Find(v2);
1392         }
1393         else {
1394           aSt2 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy2), Precision::Confusion(), N2);
1395           aStatMap.Bind(v2,aSt2);
1396           aNorMap.Bind(v2,N2);
1397         }
1398
1399         if(aNorMap.IsBound(v3)) {
1400           aSt3 = aStatMap.Find(v3);
1401           N3 = aNorMap.Find(v3);
1402         }
1403         else {
1404           aSt3 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy3), Precision::Confusion(), N3);
1405           aStatMap.Bind(v3,aSt3);
1406           aNorMap.Bind(v3,N3.XYZ());
1407         }
1408
1409         Standard_Real anAngle1 = N2.Angle(N1);
1410         Standard_Real anAngle2 = N3.Angle(N2);
1411         Standard_Real anAngle3 = N1.Angle(N3);
1412         if(aSt1 < 1 && aSt2 < 1 && aSt3 < 1 &&
1413             (anAngle1 > myAngle || anAngle2 > myAngle || anAngle3 > myAngle)) {
1414
1415           if (theIsFirst) {
1416             maxdef = -1;
1417             break;
1418           }
1419           // Record new vertex
1420           aNbPnt++;
1421           myNbLocat++;
1422           mi2d = (xy1+xy2+xy3)/3.;
1423           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1424           myLocation3d.Bind(myNbLocat,pDef);
1425           mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1426           InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1427           theInternalV.Append(InsVertex);
1428           theBadTriangles.Append(TriId);
1429         }
1430
1431         //In case if triangle is considerd as OK, we have to check three intermediate
1432         //points to be sure that we free from wave effect. If it is OK triangle passed if not split in middle point
1433         if(theBadTriangles.IsEmpty() || theBadTriangles.Last() != TriId)
1434         {
1435
1436           Bnd_Box2d aB2d;
1437           aB2d.Add(gp_Pnt2d(xy1)); aB2d.Add(gp_Pnt2d(xy2)); aB2d.Add(gp_Pnt2d(xy3));
1438
1439           Standard_Real aXMin1, aXMax1, aYMin1, aYMax1;
1440           aB2d.Get(aXMin1, aYMin1, aXMax1, aYMax1);
1441
1442           if(aXMax1 - aXMin1 < ddu && aYMax1 - aYMin1 < ddv)
1443             continue;
1444
1445           mi2d = (xy1+xy2+xy3)/3.;
1446           gp_Pnt2d aP[3] = {mi2d+(xy1-mi2d)*2/3., mi2d+(xy2-mi2d)*2/3, mi2d+(xy3-mi2d)*2/3.};
1447           gp_Dir midDir(0,0,1);
1448           Standard_Integer aSt[4];
1449           aSt[0] = GeomLib::NormEstim(BSpl, gp_Pnt2d(mi2d), Precision::Confusion(), midDir);
1450           Standard_Integer i = 0;
1451           gp_Dir dir[3] = {gp_Dir(0,0,1), gp_Dir(0,0,1), gp_Dir(0,0,1)};
1452           Standard_Real anAngle[3];
1453           for(; i < 3; i++)
1454           {
1455             gp_Dir dir(0,0,1);
1456             aSt[i+1] = GeomLib::NormEstim(BSpl, aP[i], Precision::Confusion(), dir);
1457             anAngle[i] = dir.Angle(midDir);
1458           }
1459
1460           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1461           aB2d.Add(gp_Pnt2d(xy1)); aB2d.Add(gp_Pnt2d(xy2)); aB2d.Add(gp_Pnt2d(xy3));
1462           if(aSt[0] < 1 && aSt[1] < 1 && aSt[2] < 1 && aSt[3] < 1 &&
1463             (anAngle[0] > myAngle || anAngle[1] > myAngle || anAngle[2] > myAngle) &&
1464             (aXMax1 - aXMin1 > ddu || aYMax1 - aYMin1 > ddv))
1465           {
1466             if (theIsFirst) {
1467               maxdef = -1;
1468               break;
1469             }
1470             aNbPnt++;
1471             myNbLocat++;
1472             theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1473             myLocation3d.Bind(myNbLocat,pDef);
1474             mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1475             InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1476             theInternalV.Append(InsVertex);
1477             theBadTriangles.Append(TriId);
1478           }
1479         }
1480       }
1481     }
1482
1483     if (!theIsFirst && theInternalV.Extent() > 0)
1484     {
1485       BRepMesh_Array1OfVertexOfDelaun verttab(1, theInternalV.Extent());
1486       BRepMesh_ListIteratorOfListOfVertex itVer(theInternalV);
1487       Standard_Integer ipn = 1;
1488       for (; itVer.More(); itVer.Next())
1489         verttab(ipn++) = itVer.Value();
1490
1491       theTrigu.AddVertices(verttab);
1492       nbInserted++;
1493     }
1494   }
1495   if (maxdef < 0)
1496     return maxdef;
1497   return Sqrt(maxdef);
1498 }
1499
1500 //=======================================================================
1501 //function : AddInShape
1502 //purpose  :
1503 //=======================================================================
1504 void BRepMesh_FastDiscretFace::AddInShape(const TopoDS_Face&  theFace,
1505                                           const Standard_Real theDefFace)
1506 {
1507 //  gp_Pnt Pt;
1508   BRep_Builder B;
1509   TopLoc_Location loc = theFace.Location();
1510   Handle(Poly_Triangulation) TOld = BRep_Tool::Triangulation(theFace, loc);
1511   Handle(Poly_Triangulation) TNull;
1512   Handle(Poly_PolygonOnTriangulation) NullPoly;
1513   B.UpdateFace(theFace,TNull);
1514
1515   try{
1516   BRepMesh_MapOfInteger::Iterator it;
1517
1518   Standard_Integer e1, e2, e3, nTri;
1519   Standard_Integer v1, v2, v3, iv1, iv2, iv3;
1520   Standard_Integer i, index;
1521   Standard_Boolean o1, o2, o3;
1522   TopAbs_Orientation orFace = theFace.Orientation();
1523
1524   const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1525   it.Initialize(TriMap);
1526
1527   nTri = TriMap.Extent();
1528
1529   if (nTri != 0) {
1530
1531     Poly_Array1OfTriangle Tri(1, nTri);
1532
1533     i = 1;
1534
1535     for (; it.More(); it.Next()) {
1536       myStructure->GetElement(it.Key()).Edges(e1, e2, e3, o1, o2, o3);
1537
1538       const BRepMesh_Edge& ve1=myStructure->GetLink(e1);
1539       const BRepMesh_Edge& ve2=myStructure->GetLink(e2);
1540       const BRepMesh_Edge& ve3=myStructure->GetLink(e3);
1541
1542       if (o1) {
1543         v1=ve1.FirstNode();
1544       }
1545       else {
1546         v1=ve1.LastNode();
1547       }
1548       if (o2)
1549       {
1550         v2=ve2.FirstNode();
1551         v3=ve2.LastNode();
1552       }
1553       else
1554       {
1555         v3=ve2.FirstNode();
1556         v2=ve2.LastNode();
1557       }
1558
1559       iv1 = myVemap.FindIndex(v1);
1560       if (iv1 == 0) iv1 = myVemap.Add(v1);
1561       iv2 = myVemap.FindIndex(v2);
1562       if (iv2 == 0) iv2 = myVemap.Add(v2);
1563       iv3 = myVemap.FindIndex(v3);
1564       if (iv3 == 0) iv3 = myVemap.Add(v3);
1565
1566       if (orFace == TopAbs_REVERSED) Tri(i++).Set(iv1, iv3, iv2);
1567       else Tri(i++).Set(iv1, iv2, iv3);
1568     }
1569
1570     Standard_Integer nbVertices = myVemap.Extent();
1571     Handle(Poly_Triangulation) T = new Poly_Triangulation(nbVertices, nTri, Standard_True);
1572     Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
1573     Trian = Tri;
1574     TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
1575     TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
1576
1577     for (i = 1; i <= nbVertices; i++) {
1578       index = myVemap.FindKey(i);
1579       Nodes(i) = Pnt(index);
1580       Nodes2d(i).SetXY(Vertex(index).Coord());
1581     }
1582
1583     T->Deflection(theDefFace);
1584
1585     // stockage de la triangulation dans la BRep.
1586     BRep_Builder B1;
1587     //TopLoc_Location loc = theFace.Location();
1588     if (!loc.IsIdentity()) {
1589       gp_Trsf tr = loc.Transformation();
1590       tr.Invert();
1591       for (i = Nodes.Lower(); i <= Nodes.Upper(); i++)
1592         Nodes(i).Transform(tr);
1593     }
1594     B1.UpdateFace(theFace, T);
1595
1596     // mise en place des polygones sur triangulation dans la face:
1597     BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(myInternaledges);
1598
1599     for (; It.More(); It.Next()) {
1600       const BRepMesh_PairOfPolygon& pair = It.Value();
1601       const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
1602       const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
1603       if ( NOD1 == NOD2 ) {
1604         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly, TOld,loc);
1605         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
1606       }
1607       else {
1608         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly,   TOld,loc);
1609         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
1610       }
1611     }
1612   }
1613  }
1614  catch(Standard_Failure)
1615  {
1616    //   MESH_FAILURE(theFace);
1617  }
1618 }
1619
1620
1621 //=======================================================================
1622 //function : Triangle
1623 //purpose  :
1624 //=======================================================================
1625
1626 const BRepMesh_Triangle& BRepMesh_FastDiscretFace::Triangle(const Standard_Integer Index) const
1627 {
1628   return myStructure->GetElement(Index);
1629 }
1630
1631 //=======================================================================
1632 //function : NbEdges
1633 //purpose  :
1634 //=======================================================================
1635
1636 /*Standard_Integer BRepMesh_FastDiscretFace::NbEdges() const
1637 {
1638   return myStructure->NbLinks();
1639 }*/
1640
1641 //=======================================================================
1642 //function : Edge
1643 //purpose  :
1644 //=======================================================================
1645
1646 const BRepMesh_Edge& BRepMesh_FastDiscretFace::Edge(const Standard_Integer Index) const
1647 {
1648   return myStructure->GetLink(Index);
1649 }
1650
1651
1652 //=======================================================================
1653 //function : Vertex
1654 //purpose  :
1655 //=======================================================================
1656
1657 const BRepMesh_Vertex& BRepMesh_FastDiscretFace::Vertex(const Standard_Integer Index) const
1658 {
1659   return myStructure->GetNode(Index);
1660 }
1661
1662 //=======================================================================
1663 //function : Pnt
1664 //purpose  :
1665 //=======================================================================
1666
1667 const gp_Pnt& BRepMesh_FastDiscretFace::Pnt(const Standard_Integer Index) const
1668 {
1669   return myLocation3d(myStructure->GetNode(Index).Location3d());
1670 }
1671
1672 //=======================================================================
1673 //function : FindUV
1674 //purpose  :
1675 //=======================================================================
1676
1677 gp_XY BRepMesh_FastDiscretFace::FindUV(const TopoDS_Vertex&                 theV,
1678                                        const gp_Pnt2d&                      theXY,
1679                                        const Standard_Integer               theIp,
1680                                        const Handle(BRepAdaptor_HSurface)&  theSFace,
1681                                        const Standard_Real                  theMinDist,
1682                                        BRepMesh_DataMapOfIntegerListOfXY&   theLocation2dMap)
1683 {
1684   gp_XY anUV;
1685   if (theLocation2dMap.IsBound(theIp))
1686   {
1687     BRepMesh_ListOfXY& L = theLocation2dMap.ChangeFind(theIp);
1688     anUV = L.First();
1689     if (L.Extent() != 1)
1690     {
1691       BRepMesh_ListIteratorOfListOfXY it(L);
1692       it.Next();
1693       Standard_Real dd, dmin = theXY.Distance(gp_Pnt2d(anUV));
1694       for (; it.More(); it.Next())
1695       {
1696         dd = theXY.Distance(gp_Pnt2d(it.Value()));
1697         if (dd < dmin)
1698         {
1699           anUV = it.Value();
1700           dmin = dd;
1701         }
1702       }
1703     }
1704
1705     const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(theV), theMinDist);
1706
1707     const Standard_Real Utol2d = .5 * (theSFace->LastUParameter() - theSFace->FirstUParameter());
1708     const Standard_Real Vtol2d = .5 * (theSFace->LastVParameter() - theSFace->FirstVParameter());
1709
1710     const gp_Pnt p1 = theSFace->Value(anUV.X(), anUV.Y());
1711     const gp_Pnt p2 = theSFace->Value(theXY.X(), theXY.Y());
1712
1713     if (Abs(anUV.X() - theXY.X()) > Utol2d ||
1714         Abs(anUV.Y() - theXY.Y()) > Vtol2d ||
1715         !p1.IsEqual(p2, tol))
1716     {
1717       anUV = theXY.Coord();
1718       L.Append(anUV);
1719     }
1720   }
1721   else
1722   {
1723     anUV = theXY.Coord();
1724     BRepMesh_ListOfXY L;
1725     L.Append(anUV);
1726     theLocation2dMap.Bind(theIp, L);
1727   }
1728   return anUV;
1729 }
1730
1731
1732 static Standard_Boolean GetVertexParameters(const TopoDS_Vertex& theVert,
1733                                             const TopoDS_Face& theFace,
1734                                             gp_Pnt2d& thePoint)
1735 {
1736   TopLoc_Location L;
1737   const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
1738   L = L.Predivided(theVert.Location());
1739   BRep_ListIteratorOfListOfPointRepresentation itpr =
1740     ((*((Handle(BRep_TVertex)*) &theVert.TShape()))->Points());
1741   // On regarde dabord si il y des PointRepresentation (cas non Manifold)
1742
1743   while (itpr.More()) {
1744     if (itpr.Value()->IsPointOnSurface(S,L)) {
1745       thePoint.SetCoord(itpr.Value()->Parameter(),
1746                         itpr.Value()->Parameter2());
1747       return Standard_True;
1748     }
1749     itpr.Next();
1750   }
1751   return Standard_False;
1752 }
1753
1754 //=======================================================================
1755 //function : Add
1756 //purpose  : method intended to addition internal myVertices in triangulation.
1757 //=======================================================================
1758 void BRepMesh_FastDiscretFace::Add(const TopoDS_Vertex&                theVert,
1759                                    const TopoDS_Face&                  theFace,
1760                                    const Handle(BRepAdaptor_HSurface)& thegFace)
1761 {
1762   const TopAbs_Orientation anOrient = theVert.Orientation();
1763   gp_Pnt2d uvXY;
1764   if( anOrient != TopAbs_INTERNAL || !GetVertexParameters(theVert,theFace,uvXY))
1765     return;
1766   Standard_Integer indVert =0;
1767   if (myVertices.IsBound(theVert))
1768     indVert = myVertices.Find(theVert);
1769   else
1770   {
1771     myNbLocat++;
1772     myLocation3d.Bind(myNbLocat, BRep_Tool::Pnt(theVert));
1773     indVert = myNbLocat;
1774     myVertices.Bind(theVert, indVert);
1775   }
1776   Standard_Real mindist = BRep_Tool::Tolerance(theVert);
1777   // gp_Pnt2d uvXY = BRep_Tool::Parameters(theVert,theFace);
1778   gp_XY anUV = FindUV(theVert, uvXY, indVert, thegFace, mindist, myLocation2d);
1779   BRepMesh_Vertex vf(anUV, indVert, BRepMesh_Fixed);
1780   Standard_Integer ivff = myStructure->AddNode(vf);
1781   Standard_Integer isvf = myVemap.FindIndex(ivff);
1782   if (isvf == 0) isvf = myVemap.Add(ivff);
1783 }