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