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   if( theResult.Length() == 2 )
 597   {
 598     Standard_Real    dist  = theResult.Last() - theResult.First();
 599     Standard_Integer nbint = (Standard_Integer)((dist / theFilterDist) + 0.5);
 600
 601     if( nbint > 1 )
 602     {
 603       //Five points more is maximum
 604       if( nbint > 5 )
 605       {
 606         nbint = 5;
 607       }
 608
 609       Standard_Integer i;
 610       Standard_Real dU = dist / nbint;
 611       for( i = 1; i < nbint; i++ )
 612       {
 613         theResult.InsertAfter(i, theResult.First()+i*dU);
 614       }
 615     }
 616   }
 617 }
 618
 619 void BRepMesh_FastDiscretFace::InternalVertices(const Handle(BRepAdaptor_HSurface)& theCaro,
 620                                                 BRepMesh_ListOfVertex&              theInternalV,
 621                                                 const Standard_Real                 theDefFace,
 622                                                 const BRepMesh_ClassifierPtr&       theClassifier)
 623 {
 624   BRepMesh_Vertex newV;
 625   gp_Pnt2d p2d;
 626   gp_Pnt p3d;
 627
 628   // work following the type of surface
 629   const BRepAdaptor_Surface& BS = *(BRepAdaptor_Surface*)&(theCaro->Surface());
 630   GeomAbs_SurfaceType thetype = theCaro->GetType();
 631
 632   Standard_Real umax = myAttrib->GetUMax();
 633   Standard_Real umin = myAttrib->GetUMin();
 634   Standard_Real vmax = myAttrib->GetVMax();
 635   Standard_Real vmin = myAttrib->GetVMin();
 636   Standard_Real deltaX = myAttrib->GetDeltaX();
 637   Standard_Real deltaY = myAttrib->GetDeltaY();
 638
 639   if (thetype == GeomAbs_Sphere)
 640   {
 641     gp_Sphere S = BS.Sphere();
 642     const Standard_Real R = S.Radius();
 643
 644     // Calculate parameters for iteration in V direction
 645     Standard_Real Dv = 1.0 - (theDefFace/R);
 646     if (Dv < 0.0) Dv = 0.0;
 647     Standard_Real oldDv = 2.0 * ACos (Dv);
 648     Dv  =  .7 * oldDv; //.7 ~= sqrt(2.) - Dv is hypotenuse of triangle when oldDv is legs
 649     const Standard_Real sv = vmax - vmin;
 650     Dv = sv/((Standard_Integer)(sv/Dv) + 1);
 651     const Standard_Real pasvmax = vmax-Dv*0.5;
 652
 653     //Du can be defined from relation: 2*r*Sin(Du/2) = 2*R*Sin(Dv/2), r = R*Cos(v)
 654     //here approximate relation r*Du = R*Dv is used
 655
 656     Standard_Real Du, pasu, pasv; //, ru;
 657     const Standard_Real su = umax-umin;
 658     Standard_Boolean Shift = Standard_False;
 659     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv)
 660     {
 661       // Calculate parameters for iteration in U direction
 662       // 1.-.365*pasv*pasv is simple approximation of Cos(pasv)
 663       // with condition that it gives ~.1 when pasv = pi/2
 664       Du = Dv/(1.-.365*pasv*pasv);
 665       Du = su/((Standard_Integer)(su/Du) + 1);
 666       Shift = !Shift;
 667       const Standard_Real d = (Shift)? Du*.5 : 0.;
 668       const Standard_Real pasumax = umax-Du*0.5 + d;
 669       for (pasu = umin + Du - d; pasu < pasumax; pasu += Du)
 670       {
 671         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 672         {
 673           // Record 3d point
 674 #ifdef DEB_MESH_CHRONO
 675           D0Internal++;
 676 #endif
 677           ElSLib::D0(pasu, pasv, S, p3d);
 678           myNbLocat++;
 679           myLocation3d.Bind(myNbLocat, p3d);
 680           // Record 2d point
 681           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 682           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 683           theInternalV.Append(newV);
 684         }
 685       }
 686     }
 687   }
 688   else if (thetype == GeomAbs_Cylinder)
 689   {
 690     gp_Cylinder S = BS.Cylinder();
 691     const Standard_Real R = S.Radius();
 692
 693     // Calculate parameters for iteration in U direction
 694     Standard_Real Du = 1.0 - (theDefFace/R);
 695     if (Du < 0.0) Du = 0.0;
 696     Du = 2.0 * ACos (Du);
 697     if (Du > myAngle) Du = myAngle;
 698     const Standard_Real su = umax - umin;
 699     const Standard_Integer nbU = (Standard_Integer)(su/Du);
 700     Du = su/(nbU+1);
 701
 702     // Calculate parameters for iteration in V direction
 703     const Standard_Real sv = vmax - vmin;
 704     Standard_Integer nbV = (Standard_Integer)( nbU*sv/(su*R) );
 705     nbV = Min(nbV, 100*nbU);
 706     Standard_Real Dv = sv/(nbV+1);
 707
 708     Standard_Real pasu, pasv, pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
 709     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
 710       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
 711         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 712         {
 713           // Record 3d point
 714           ElSLib::D0(pasu, pasv, S, p3d);
 715           myNbLocat++;
 716           myLocation3d.Bind(myNbLocat, p3d);
 717           // Record 2d point
 718           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 719           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 720           theInternalV.Append(newV);
 721         }
 722       }
 723     }
 724   }
 725   else if (thetype == GeomAbs_Cone)
 726   {
 727     Standard_Real R, RefR, SAng;
 728     gp_Cone C = BS.Cone();
 729     RefR = C.RefRadius();
 730     SAng = C.SemiAngle();
 731     R = Max(Abs(RefR+vmin*Sin(SAng)), Abs(RefR+vmax*Sin(SAng)));
 732     Standard_Real Du, Dv, pasu, pasv;
 733     Du = Max(1.0e0 - (theDefFace/R),0.0e0);
 734     Du  = (2.0 * ACos (Du));
 735     Standard_Integer nbU = (Standard_Integer) ( (umax-umin)/Du );
 736     Standard_Integer nbV = (Standard_Integer) ( nbU*(vmax-vmin)/((umax-umin)*R) );
 737     Du = (umax-umin)/(nbU+1);
 738     Dv = (vmax-vmin)/(nbV+1);
 739
 740     Standard_Real pasvmax = vmax-Dv*0.5, pasumax = umax-Du*0.5;
 741     for (pasv = vmin + Dv; pasv < pasvmax; pasv += Dv) {
 742       for (pasu = umin + Du; pasu < pasumax; pasu += Du) {
 743         if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 744         {
 745           // Record 3d point
 746           ElSLib::D0(pasu, pasv, C, p3d);
 747           myNbLocat++;
 748           myLocation3d.Bind(myNbLocat, p3d);
 749           // Record 2d point
 750           p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 751           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 752           theInternalV.Append(newV);
 753         }
 754       }
 755     }
 756   }
 757   else if (thetype == GeomAbs_Torus)
 758   {
 759     gp_Torus T = BS.Torus();
 760
 761     Standard_Boolean insert;
 762     Standard_Integer i, j, ParamULength, ParamVLength;
 763     Standard_Real pp, pasu, pasv;
 764     Standard_Real r = T.MinorRadius(), R = T.MajorRadius();
 765
 766     TColStd_SequenceOfReal ParamU, ParamV;
 767
 768     Standard_Real Du, Dv;//, pasu, pasv;
 769     Dv = Max(1.0e0 - (theDefFace/r),0.0e0) ;
 770     Standard_Real oldDv = 2.0 * ACos (Dv);
 771     oldDv = Min(oldDv, myAngle);
 772     Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
 773     Dv = oldDv;
 774
 775     Standard_Integer nbV = Max((Standard_Integer)((vmax-vmin)/Dv), 2);
 776     Dv = (vmax-vmin)/(nbV+1);
 777     Standard_Real ru = R + r;
 778     if (ru > 1.e-16)
 779     {
 780       Du  = 2.0 * ACos(Max(1.0 - (theDefFace/ru),0.0));
 781       if (myAngle < Du) Du = myAngle;
 782       Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
 783       if(aa < gp::Resolution())
 784         return;
 785       Du *= Min(oldDv, Du) / aa;
 786     }
 787     else Du = Dv;
 788
 789     Standard_Integer nbU = Max((Standard_Integer)((umax-umin)/Du), 2);
 790     nbU = Max(nbU , (int)(nbV*(umax-umin)*R/((vmax-vmin)*r)/5.));
 791     Du = (umax-umin)/(nbU+1);
 792
 793     if (R < r)
 794     {
 795       // As the points of edges are returned.
 796       // in this case, the points are not representative.
 797
 798       //-- Choose DeltaX and DeltaY so that to avoid skipping points on the grid
 799            for (i = 0; i <= nbU; i++) ParamU.Append(umin + i* Du);
 800     }//R<r
 801     else //U if R > r
 802     {
 803       //--ofv: U
 804       // Number of mapped U parameters
 805       const Standard_Integer LenU = myUParam.Extent();
 806       // Fill array of U parameters
 807       TColStd_Array1OfReal Up(1,LenU);
 808       for (j = 1; j <= LenU; j++) Up(j) = myUParam(j);
 809
 810       // Calculate DU, leave array of parameters
 811       Standard_Real aDU = FUN_CalcAverageDUV(Up,LenU);
 812       aDU = Max(aDU, Abs(umax -  umin) / (Standard_Real) nbU / 2.);
 813       Standard_Real dUstd = Abs(umax -  umin) / (Standard_Real) LenU;
 814       if (aDU > dUstd) dUstd = aDU;
 815       // Add U parameters
 816       for (j = 1; j <= LenU; j++)
 817       {
 818         pp = Up(j);
 819         insert = Standard_True;
 820         ParamULength = ParamU.Length();
 821         for (i = 1; i <= ParamULength && insert; i++)
 822         {
 823           insert = (Abs(ParamU.Value(i)-pp) > (0.5*dUstd));
 824         }
 825         if (insert) ParamU.Append(pp);
 826       }
 827     }
 828
 829     //--ofv: V
 830     // Number of mapped V parameters
 831     const Standard_Integer LenV = myVParam.Extent();
 832     // Fill array of V parameters
 833     TColStd_Array1OfReal Vp(1,LenV);
 834     for (j = 1; j <= LenV; j++) Vp(j) = myVParam(j);
 835     // Calculate DV, sort array of parameters
 836     Standard_Real aDV = FUN_CalcAverageDUV(Vp,LenV);
 837     aDV = Max(aDV, Abs(vmax -  vmin) / (Standard_Real) nbV / 2.);
 838
 839     Standard_Real dVstd = Abs(vmax -  vmin) / (Standard_Real) LenV;
 840     if (aDV > dVstd) dVstd = aDV;
 841     // Add V parameters
 842     for (j = 1; j <= LenV; j++)
 843     {
 844       pp = Vp(j);
 845
 846       insert = Standard_True;
 847       ParamVLength = ParamV.Length();
 848       for (i = 1; i <= ParamVLength && insert; i++)
 849       {
 850         insert = (Abs(ParamV.Value(i)-pp) > (dVstd*2./3.));
 851       }
 852       if (insert) ParamV.Append(pp);
 853     }
 854
 855     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
 856     Standard_Real uminnew = umin+deltaY*0.1;
 857     Standard_Real vminnew = vmin+deltaX*0.1;
 858     Standard_Real umaxnew = umax-deltaY*0.1;
 859     Standard_Real vmaxnew = vmax-deltaX*0.1;
 860
 861     for (i = 1; i <= Lu; i++)
 862     {
 863       pasu = ParamU.Value(i);
 864       if (pasu >= uminnew && pasu < umaxnew)
 865       {
 866         for (j = 1; j <= Lv; j++)
 867         {
 868           pasv = ParamV.Value(j);
 869           if (pasv >= vminnew && pasv < vmaxnew)
 870           {
 871             if (theClassifier->Perform(gp_Pnt2d(pasu, pasv)) == TopAbs_IN)
 872             {
 873               // Record 3d point
 874               ElSLib::D0(pasu, pasv, T, p3d);
 875               myNbLocat++;
 876               myLocation3d.Bind(myNbLocat, p3d);
 877               // Record 2d point
 878               p2d.SetCoord((pasu-umin)/deltaX, (pasv-vmin)/deltaY);
 879               newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
 880               theInternalV.Append(newV);
 881             }
 882           }
 883         }
 884       }
 885     }
 886   }
 887   else if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
 888   {
 889     //define resolutions
 890     Standard_Real uRes = BS.UResolution(theDefFace);
 891     Standard_Real vRes = BS.VResolution(theDefFace);
 892
 893     // Sort and filter sequence of parameters
 894     Standard_Real aMinDu = Precision::PConfusion();
 895     if(deltaX < 1.)
 896       aMinDu/=deltaX;
 897
 898     Standard_Real aDuMaxLim = 0.1*(umax-umin);
 899     Standard_Real ddu = Min(aDuMaxLim,Max(0.005*(umax-umin),2.*uRes));
 900     TColStd_SequenceOfReal ParamU;
 901     filterParameters(myUParam,aMinDu,ddu,ParamU);
 902     Standard_Integer ParamULength = ParamU.Length();
 903
 904     Standard_Real aMinDv = Precision::PConfusion();
 905     if(deltaY < 1)
 906       aMinDv/=deltaY;
 907
 908     Standard_Real aDvMaxLim = 0.1*(vmax-vmin);
 909     Standard_Real ddv = Min(aDvMaxLim,Max(0.005*(vmax-vmin),2.*vRes));
 910     TColStd_SequenceOfReal ParamV;
 911     filterParameters(myVParam,aMinDv,ddv,ParamV);
 912     Standard_Integer ParamVLength = ParamV.Length();
 913
 914     // check intermediate isolines
 915     Handle(Geom_Surface) B;
 916     if (thetype == GeomAbs_BezierSurface) {
 917       B = BS.Bezier();
 918     }
 919     else {
 920       B = BS.BSpline();
 921     }
 922
 923     gp_Pnt P1, P2, PControl;
 924     Standard_Real u, v, dist;
 925
 926     // precision for compare square distances
 927     double dPreci = Precision::Confusion()*Precision::Confusion();
 928
 929     // Insert V parameters by deflection criterion
 930     Standard_Integer i,j;
 931     Standard_Real V1, V2, U1, U2;
 932     for (i = 1; i <= ParamULength; i++) {
 933       Handle(Geom_Curve) IsoU = B->UIso(ParamU.Value(i));
 934       V1 = ParamV.Value(1);
 935       P1 = IsoU->Value(V1);
 936       for (j = 2; j <= ParamVLength;) {
 937         V2 = ParamV.Value(j);
 938         P2 = IsoU->Value(V2);
 939         v = 0.5*(V1+V2);
 940         PControl = IsoU->Value(v);
 941         // 23.03.2010 skl for OCC21645 - change precision for comparison
 942         if( P1.SquareDistance(P2) > dPreci ) {
 943           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
 944           dist = L.Distance(PControl);
 945         }
 946         else {
 947           dist = P1.Distance(PControl);
 948         }
 949         if (dist > theDefFace) {
 950           // insertion
 951           ParamV.InsertBefore(j, v);
 952           ParamVLength++;
 953         }
 954         else {
 955           //put regular grig for normals
 956           gp_Dir N1(0,0,1),N2(0,0,1);
 957           Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),V1), Precision::Confusion(), N1);
 958           Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(ParamU.Value(i),v), Precision::Confusion(), N2);
 959
 960           Standard_Real anAngle1 = N2.Angle(N1);
 961           if(aSt1 < 1 && aSt2 < 1 && anAngle1 > myAngle ) {
 962             // insertion
 963             ParamV.InsertBefore(j, v);
 964             ParamVLength++;
 965           }
 966           else {
 967             V1 = V2;
 968             P1 = P2;
 969             j++;
 970           }
 971         }
 972       }
 973     }
 974
 975     for (i = 2; i < ParamVLength; i++) {
 976       v = ParamV.Value(i);
 977       Handle(Geom_Curve) IsoV = B->VIso(v);
 978       U1 = ParamU.Value(1);
 979       P1 = IsoV->Value(U1);
 980       for (j = 2; j <= ParamULength;) {
 981         U2 = ParamU.Value(j);
 982         P2 = IsoV->Value(U2);
 983         u = 0.5*(U1+U2);
 984         PControl = IsoV->Value(u);
 985         // 23.03.2010 skl for OCC21645 - change precision for comparison
 986         if( P1.SquareDistance(P2) > dPreci ) {
 987           gp_Lin L (P1, gp_Dir(gp_Vec(P1, P2)));
 988           dist = L.Distance(PControl);
 989         }
 990         else {
 991           dist = P1.Distance(PControl);
 992         }
 993         if (dist > theDefFace) {
 994           // insertion
 995           ParamU.InsertBefore(j, u);
 996           ParamULength++;
 997         }
 998         else {
 999           //check normal
1000           //put regular grig for normals
1001           gp_Dir N1(0,0,1),N2(0,0,1);
1002           Standard_Boolean aSt1 = GeomLib::NormEstim(B, gp_Pnt2d(U1,v), Precision::Confusion(), N1);
1003           Standard_Boolean aSt2 = GeomLib::NormEstim(B, gp_Pnt2d(u,v), Precision::Confusion(), N2);
1004
1005           Standard_Real anAngle1 = N2.Angle(N1);
1006                 if(aSt1 < 1 && aSt2 < 1 && anAngle1 > myAngle) {
1007             // insertion
1008             ParamU.InsertBefore(j, u);
1009             ParamULength++;
1010           }
1011           else {
1012             U1 = U2;
1013             P1 = P2;
1014             if (j < ParamULength) {
1015               // Classify intersection point
1016               if (theClassifier->Perform(gp_Pnt2d(U1, v)) == TopAbs_IN)
1017               {
1018                 // Record 3d point
1019                 myNbLocat++;
1020                 myLocation3d.Bind(myNbLocat, P1);
1021                 // Record 2d point
1022                 p2d.SetCoord((U1-umin)/deltaX, (v-vmin)/deltaY);
1023                 newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
1024                 theInternalV.Append(newV);
1025               }
1026             }
1027             j++;
1028           }
1029         }
1030       }
1031     }
1032   }
1033   else {
1034     const Standard_Real anAngle = 0.35;
1035
1036     Standard_Integer i, j, nbpointsU = 10, nbpointsV = 10;
1037     Adaptor3d_IsoCurve tabu[11], tabv[11];
1038
1039     TColStd_SequenceOfReal ParamU, ParamV;
1040     Standard_Real u, v, du, dv;
1041     Standard_Integer iu, iv;
1042     Standard_Real f, l;
1043
1044     du = (umax-umin) / (nbpointsU+1); dv = (vmax-vmin) / (nbpointsV+1);
1045
1046     for (iu = 0; iu <= nbpointsU; iu++) {
1047       u = umin + iu*du;
1048       tabu[iu].Load(theCaro);
1049       tabu[iu].Load(GeomAbs_IsoU, u);
1050     }
1051
1052     for (iv = 0; iv <= nbpointsV; iv++) {
1053       v = vmin + iv*dv;
1054       tabv[iv].Load(theCaro);
1055       tabv[iv].Load(GeomAbs_IsoV, v);
1056     }
1057
1058     Standard_Integer imax = 1, MaxV = 0;
1059
1060     GCPnts_TangentialDeflection* tabGU = new GCPnts_TangentialDeflection[nbpointsU+1];
1061
1062     for (i = 0; i <= nbpointsU; i++) {
1063       f = Max(vmin, tabu[i].FirstParameter());
1064       l = Min(vmax, tabu[i].LastParameter());
1065       GCPnts_TangentialDeflection theDeflection(tabu[i], f, l, anAngle, 0.7*theDefFace, 2);
1066       tabGU[i] = theDeflection;
1067       if (tabGU[i].NbPoints() > MaxV) {
1068         MaxV = tabGU[i].NbPoints();
1069         imax = i;
1070       }
1071     }
1072
1073     // return table of parameters V:
1074     Standard_Integer NV = tabGU[imax].NbPoints();
1075     for (i = 1; i <= NV; i++) {
1076       ParamV.Append(tabGU[imax].Parameter(i));
1077     }
1078     delete [] tabGU;
1079
1080     imax = 1;
1081     Standard_Integer MaxU = 0;
1082
1083     GCPnts_TangentialDeflection* tabGV = new GCPnts_TangentialDeflection[nbpointsV+1];
1084
1085     for (i = 0; i <= nbpointsV; i++) {
1086       f = Max(umin, tabv[i].FirstParameter());
1087       l = Min(umax, tabv[i].LastParameter());
1088       GCPnts_TangentialDeflection thedeflection2(tabv[i], f, l, anAngle, 0.7*theDefFace, 2);
1089       tabGV[i] = thedeflection2;
1090       if (tabGV[i].NbPoints() > MaxU) {
1091         MaxU = tabGV[i].NbPoints();
1092         imax = i;
1093       }
1094     }
1095
1096     // return table of parameters U:
1097     Standard_Integer NU = tabGV[imax].NbPoints();
1098     for (i = 1; i <= NU; i++) {
1099       ParamU.Append(tabGV[imax].Parameter(i));
1100     }
1101     delete [] tabGV;
1102
1103     if (ParamU.Length() == 2) {
1104       ParamU.InsertAfter(1, (umax+umin)*0.5);
1105     }
1106     if (ParamV.Length() == 2) {
1107       ParamV.InsertAfter(1, (vmax+vmin)*0.5);
1108     }
1109
1110     TColStd_SequenceOfReal InsertV, InsertU;
1111     gp_Pnt P1;
1112
1113     Adaptor3d_IsoCurve IsoV;
1114     IsoV.Load(theCaro);
1115
1116     Standard_Integer Lu = ParamU.Length(), Lv = ParamV.Length();
1117
1118     for (i = 2; i < Lv; i++) {
1119       v = ParamV.Value(i);
1120       IsoV.Load(GeomAbs_IsoV, v);
1121       for (j = 2; j < Lu; j++) {
1122         u = ParamU.Value(j);
1123         if (theClassifier->Perform(gp_Pnt2d(u, v)) == TopAbs_IN)
1124         {
1125           // Record 3d point
1126           P1 = IsoV.Value(u);
1127           myNbLocat++;
1128           myLocation3d.Bind(myNbLocat, P1);
1129           // Record 2d point
1130           p2d.SetCoord((u-umin)/deltaX, (v-vmin)/deltaY);
1131           newV.Initialize(p2d.XY(), myNbLocat, BRepMesh_Free);
1132           theInternalV.Append(newV);
1133         }
1134       }
1135     }
1136   }
1137 #ifdef DEB_MESH_CHRONO
1138   chInternal.Stop();
1139 #endif
1140 }
1141
1142 /**
1143 *  Internal class Couple, moved from MeshData package
1144 */
1145
1146 class BRepMesh_Couple
1147 {
1148  public:
1149   BRepMesh_Couple() { myI1 = myI2 = 0; }
1150   BRepMesh_Couple(const Standard_Integer I1,
1151                   const Standard_Integer I2)
1152   { myI1 = I1; myI2 = I2; }
1153
1154   Standard_Integer myI1;
1155   Standard_Integer myI2;
1156 };
1157
1158 inline Standard_Boolean IsEqual(const BRepMesh_Couple& one,
1159                                 const BRepMesh_Couple& other)
1160 {
1161   if (one.myI1 == other.myI1 &&
1162       one.myI2 == other.myI2) return Standard_True;
1163   else return Standard_False;
1164 }
1165
1166 inline Standard_Integer HashCode(const BRepMesh_Couple& one,
1167                                  const Standard_Integer Upper)
1168 {
1169   return ::HashCode((one.myI1+one.myI2), Upper);
1170 }
1171
1172 typedef NCollection_Map<BRepMesh_Couple> BRepMesh_MapOfCouple;
1173
1174 //=======================================================================
1175 //function : Control
1176 //purpose  :
1177 //=======================================================================
1178 Standard_Real BRepMesh_FastDiscretFace::Control(const Handle(BRepAdaptor_HSurface)& theCaro,
1179                                                 const Standard_Real                 theDefFace,
1180                                                 BRepMesh_ListOfVertex&              theInternalV,
1181                                                 TColStd_ListOfInteger&              theBadTriangles,
1182                                                 TColStd_ListOfInteger&              theNulTriangles,
1183                                                 BRepMesh_Delaun&                    theTrigu,
1184                                                 const Standard_Boolean              theIsFirst)
1185 {
1186   //IMPORTANT: Constants used in calculations
1187   const Standard_Real MinimalArea2d = 1.e-9;
1188   const Standard_Real MinimalSqLength3d = 1.e-12;
1189
1190   // Define the number of iterations
1191   Standard_Integer myNbIterations = 11;
1192   const Standard_Integer nbPasses = (theIsFirst? 1 : myNbIterations);
1193
1194   // Initialize stop condition
1195   Standard_Boolean allDegenerated = Standard_False;
1196   Standard_Integer nbInserted = 1;
1197
1198   // Create map of links to skip already processed
1199   Standard_Integer nbtriangles;
1200
1201   nbtriangles = myStructure->ElemOfDomain().Extent();
1202   if (nbtriangles <= 0) return -1.0;
1203   BRepMesh_MapOfCouple theCouples(3*nbtriangles);
1204
1205   gp_XY mi2d;
1206   gp_XYZ vecEd1, vecEd2, vecEd3;
1207   gp_Pnt pDef;
1208   Standard_Real dv = 0, defl = 0, maxdef = -1;
1209   Standard_Integer pass = 1, nf = 0, nl = 0;
1210   BRepMesh_Vertex InsVertex;
1211   Standard_Boolean caninsert;
1212
1213   Standard_Real sqdefface = theDefFace * theDefFace;
1214
1215   GeomAbs_SurfaceType thetype = theCaro->GetType();
1216   Handle(Geom_Surface) BSpl;
1217   Standard_Boolean isSpline = Standard_False;
1218   if (thetype == GeomAbs_BezierSurface || thetype == GeomAbs_BSplineSurface)
1219   {
1220     isSpline = Standard_True;
1221     if (thetype == GeomAbs_BezierSurface)
1222       BSpl = theCaro->Bezier();
1223     else
1224       BSpl = theCaro->BSpline();
1225   }
1226
1227   NCollection_DataMap<Standard_Integer,gp_Dir> aNorMap;
1228   NCollection_DataMap<Standard_Integer,Standard_Integer> aStatMap;
1229
1230   // Perform refinement passes
1231   for (; pass <= nbPasses && nbInserted && !allDegenerated; pass++)
1232   {
1233     theInternalV.Clear();
1234     theBadTriangles.Clear();
1235
1236     // Reset stop condition
1237     allDegenerated = Standard_True;
1238     nbInserted = 0;
1239     maxdef = -1.0;
1240
1241     // Do not insert nodes in last pass in non-SharedMode
1242     caninsert = (myWithShare || pass < nbPasses);
1243
1244     // Read mesh size
1245     nbtriangles = myStructure->ElemOfDomain().Extent();
1246     if (nbtriangles <= 0) break;
1247
1248     // Iterate on current triangles
1249     BRepMesh_MapOfInteger::Iterator triDom;
1250     const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1251     triDom.Initialize(TriMap);
1252     Standard_Integer aNbPnt = 0;
1253     Standard_Real umin = myAttrib->GetUMin();
1254     Standard_Real vmin = myAttrib->GetVMin();
1255     Standard_Real deltaX = myAttrib->GetDeltaX();
1256     Standard_Real deltaY = myAttrib->GetDeltaY();
1257     for (; triDom.More(); triDom.Next())
1258     {
1259       Standard_Integer TriId = triDom.Key();
1260       const BRepMesh_Triangle& curTri=Triangle(TriId);
1261       if (curTri.Movability()==BRepMesh_Deleted) continue;
1262
1263       Standard_Boolean o1, o2, o3;
1264       Standard_Integer v1 = 0, v2 = 0, v3 = 0, e1 = 0, e2 = 0, e3 = 0;
1265       curTri.Edges(e1, e2, e3, o1, o2, o3);
1266
1267       const BRepMesh_Edge& edg1=Edge(e1);
1268       const BRepMesh_Edge& edg2=Edge(e2);
1269       const BRepMesh_Edge& edg3=Edge(e3);
1270
1271       Standard_Boolean m1 = (edg1.Movability() == BRepMesh_Frontier);
1272       Standard_Boolean m2 = (edg2.Movability() == BRepMesh_Frontier);
1273       Standard_Boolean m3 = (edg3.Movability() == BRepMesh_Frontier);
1274       if (o1) {
1275         v1=edg1.FirstNode();
1276         v2=edg1.LastNode();
1277       }
1278       else {
1279         v1=edg1.LastNode();
1280         v2=edg1.FirstNode();
1281       }
1282       if (o2)
1283         v3=edg2.LastNode();
1284       else
1285         v3=edg2.FirstNode();
1286
1287       const BRepMesh_Vertex& vert1=Vertex(v1);
1288       const BRepMesh_Vertex& vert2=Vertex(v2);
1289       const BRepMesh_Vertex& vert3=Vertex(v3);
1290
1291       const gp_XYZ& p1=myLocation3d(vert1.Location3d()).Coord();
1292       const gp_XYZ& p2=myLocation3d(vert2.Location3d()).Coord();
1293       const gp_XYZ& p3=myLocation3d(vert3.Location3d()).Coord();
1294
1295       vecEd1 = p2 - p1;
1296       vecEd2 = p3 - p2;
1297       vecEd3 = p1 - p3;
1298
1299       // Check for degenerated triangle
1300       if (vecEd1.SquareModulus() < MinimalSqLength3d ||
1301           vecEd2.SquareModulus() < MinimalSqLength3d ||
1302           vecEd3.SquareModulus() < MinimalSqLength3d)
1303       {
1304         theNulTriangles.Append(TriId);
1305         continue;
1306       }
1307
1308       allDegenerated = Standard_False;
1309
1310       gp_XY xy1(vert1.Coord().X()*deltaX+umin,vert1.Coord().Y()*deltaY+vmin);
1311       gp_XY xy2(vert2.Coord().X()*deltaX+umin,vert2.Coord().Y()*deltaY+vmin);
1312       gp_XY xy3(vert3.Coord().X()*deltaX+umin,vert3.Coord().Y()*deltaY+vmin);
1313
1314       // Check triangle area in 2d
1315       if (Abs((xy2-xy1)^(xy3-xy1)) < MinimalArea2d)
1316       {
1317         theNulTriangles.Append(TriId);
1318         continue;
1319       }
1320
1321       // Check triangle normal
1322       gp_XYZ normal(vecEd1^vecEd2);
1323       dv = normal.Modulus();
1324       if (dv < Precision::Confusion())
1325       {
1326         theNulTriangles.Append(TriId);
1327         continue;
1328       }
1329       normal /= dv;
1330
1331       // Check deflection on triangle
1332       mi2d = (xy1+xy2+xy3)/3.0;
1333       theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1334       defl = Abs(normal*(pDef.XYZ()-p1));
1335       defl = defl*defl;
1336       /*mi2d = (xy1+xy2+xy3)/3.0;
1337       theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1338       defl = pDef.SquareDistance((p1+p2+p3)/3.);*/
1339       if (defl > maxdef) maxdef = defl;
1340       if (defl > sqdefface)
1341       {
1342         if (theIsFirst) break;
1343         if (caninsert)
1344         {
1345           // Record new vertex
1346           aNbPnt++;
1347           myNbLocat++;
1348           myLocation3d.Bind(myNbLocat,pDef);
1349           mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1350           InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1351           theInternalV.Append(InsVertex);
1352         }
1353         theBadTriangles.Append(TriId);
1354       }
1355
1356       if (!m2) // Not a boundary
1357       {
1358         // Check if this link was already processed
1359         if (v2 < v3) { nf = v2; nl = v3; } else { nf = v3; nl = v2; }
1360         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1361         {
1362           // Check deflection on edge 1
1363           mi2d = (xy2+xy3)*0.5;
1364           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1365           gp_Lin L (p2, gp_Vec(p2, p3));
1366           defl = L.SquareDistance(pDef);
1367           if (defl > maxdef) maxdef = defl;
1368           if (defl > sqdefface)
1369           {
1370             if (theIsFirst) break;
1371             if (caninsert)
1372             {
1373               // Record new vertex
1374               aNbPnt++;
1375               myNbLocat++;
1376               myLocation3d.Bind(myNbLocat,pDef);
1377               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1378               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1379               theInternalV.Append(InsVertex);
1380             }
1381             theBadTriangles.Append(TriId);
1382           }
1383         }
1384       }
1385
1386       if (!m3) // Not a boundary
1387       {
1388         // Check if this link was already processed
1389         if (v1 < v3) { nf = v1; nl = v3; } else { nf = v3; nl = v1; }
1390         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1391         {
1392           // Check deflection on edge 2
1393           mi2d = (xy3+xy1)*0.5;
1394           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1395           gp_Lin L (p1, gp_Vec(p1, p3));
1396           defl = L.SquareDistance(pDef);
1397           if (defl > maxdef) maxdef = defl;
1398           if (defl > sqdefface)
1399           {
1400             if (theIsFirst) break;
1401             if (caninsert)
1402             {
1403               // Record new vertex
1404               aNbPnt++;
1405               myNbLocat++;
1406               myLocation3d.Bind(myNbLocat,pDef);
1407               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1408               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1409               theInternalV.Append(InsVertex);
1410             }
1411             theBadTriangles.Append(TriId);
1412           }
1413         }
1414       }
1415
1416       if (!m1) // Not a boundary
1417       {
1418         // Check if this link was already processed
1419         if (v1 < v2) { nf = v1; nl = v2; } else { nf = v2; nl = v1; }
1420         if (theCouples.Add(BRepMesh_Couple(nf,nl)))
1421         {
1422           // Check deflection on edge 3
1423           mi2d = (xy1+xy2)*0.5;
1424           theCaro->D0(mi2d.X(), mi2d.Y(), pDef);
1425           gp_Lin L (p1, gp_Vec(p1, p2));
1426           defl = L.SquareDistance(pDef);
1427           if (defl > maxdef) maxdef = defl;
1428           if (defl > sqdefface)
1429           {
1430             if (theIsFirst) break;
1431             if (caninsert)
1432             {
1433               // Record new vertex
1434               aNbPnt++;
1435               myNbLocat++;
1436               myLocation3d.Bind(myNbLocat,pDef);
1437               mi2d.SetCoord((mi2d.X()-umin)/deltaX,(mi2d.Y()-vmin)/deltaY);
1438               InsVertex.Initialize(mi2d,myNbLocat,BRepMesh_Free);
1439               theInternalV.Append(InsVertex);
1440             }
1441             theBadTriangles.Append(TriId);
1442           }
1443         }
1444       }
1445
1446       //check normal on bsplines
1447       if(theIsFirst && isSpline && !BSpl.IsNull() )
1448       {
1449         gp_Dir N1(0,0,1), N2(0,0,1), N3(0,0,1);
1450         Standard_Integer aSt1, aSt2, aSt3;
1451         if(aNorMap.IsBound(v1)) {
1452           aSt1 = aStatMap.Find(v1);
1453           N1 =aNorMap.Find(v1);
1454         }
1455         else {
1456           aSt1 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy1), Precision::Confusion(), N1);
1457           aStatMap.Bind(v1,aSt1);
1458           aNorMap.Bind(v1,N1);
1459         }
1460
1461         if(aNorMap.IsBound(v2)) {
1462           aSt2 = aStatMap.Find(v2);
1463           N2 = aNorMap.Find(v2);
1464         }
1465         else {
1466           aSt2 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy2), Precision::Confusion(), N2);
1467           aStatMap.Bind(v2,aSt2);
1468           aNorMap.Bind(v2,N2);
1469         }
1470
1471         if(aNorMap.IsBound(v3)) {
1472           aSt3 = aStatMap.Find(v3);
1473           N3 = aNorMap.Find(v3);
1474         }
1475         else {
1476           aSt3 = GeomLib::NormEstim(BSpl, gp_Pnt2d(xy3), Precision::Confusion(), N3);
1477           aStatMap.Bind(v3,aSt3);
1478           aNorMap.Bind(v3,N3.XYZ());
1479         }
1480
1481         Standard_Real anAngle1 = N2.Angle(N1);
1482         Standard_Real anAngle2 = N3.Angle(N2);
1483         Standard_Real anAngle3 = N1.Angle(N3);
1484         if(aSt1 < 1 && aSt2 < 1 && aSt3 < 1 &&
1485             (anAngle1 > myAngle || anAngle2 > myAngle || anAngle3 > myAngle)) {
1486
1487             maxdef = -1;
1488             break;
1489         }
1490       }
1491     }
1492
1493     if (!theIsFirst && theInternalV.Extent() > 0)
1494     {
1495       BRepMesh_Array1OfVertexOfDelaun verttab(1, theInternalV.Extent());
1496       BRepMesh_ListIteratorOfListOfVertex itVer(theInternalV);
1497       Standard_Integer ipn = 1;
1498       for (; itVer.More(); itVer.Next())
1499         verttab(ipn++) = itVer.Value();
1500
1501       theTrigu.AddVertices(verttab);
1502       nbInserted++;
1503     }
1504   }
1505
1506    if (maxdef < 0)
1507     return maxdef;
1508   return Sqrt(maxdef);
1509 }
1510
1511 //=======================================================================
1512 //function : AddInShape
1513 //purpose  :
1514 //=======================================================================
1515 void BRepMesh_FastDiscretFace::AddInShape(const TopoDS_Face&  theFace,
1516                                           const Standard_Real theDefFace)
1517 {
1518 //  gp_Pnt Pt;
1519   BRep_Builder B;
1520   TopLoc_Location loc = theFace.Location();
1521   Handle(Poly_Triangulation) TOld = BRep_Tool::Triangulation(theFace, loc);
1522   Handle(Poly_Triangulation) TNull;
1523   Handle(Poly_PolygonOnTriangulation) NullPoly;
1524   B.UpdateFace(theFace,TNull);
1525
1526   try{
1527   BRepMesh_MapOfInteger::Iterator it;
1528
1529   Standard_Integer e1, e2, e3, nTri;
1530   Standard_Integer v1, v2, v3, iv1, iv2, iv3;
1531   Standard_Integer i, index;
1532   Standard_Boolean o1, o2, o3;
1533   TopAbs_Orientation orFace = theFace.Orientation();
1534
1535   const BRepMesh_MapOfInteger& TriMap = myStructure->ElemOfDomain();
1536   it.Initialize(TriMap);
1537
1538   nTri = TriMap.Extent();
1539
1540   if (nTri != 0) {
1541
1542     Poly_Array1OfTriangle Tri(1, nTri);
1543
1544     i = 1;
1545
1546     for (; it.More(); it.Next()) {
1547       myStructure->GetElement(it.Key()).Edges(e1, e2, e3, o1, o2, o3);
1548
1549       const BRepMesh_Edge& ve1=myStructure->GetLink(e1);
1550       const BRepMesh_Edge& ve2=myStructure->GetLink(e2);
1551
1552       if (o1) {
1553         v1=ve1.FirstNode();
1554       }
1555       else {
1556         v1=ve1.LastNode();
1557       }
1558       if (o2)
1559       {
1560         v2=ve2.FirstNode();
1561         v3=ve2.LastNode();
1562       }
1563       else
1564       {
1565         v3=ve2.FirstNode();
1566         v2=ve2.LastNode();
1567       }
1568
1569       iv1 = myVemap.FindIndex(v1);
1570       if (iv1 == 0) iv1 = myVemap.Add(v1);
1571       iv2 = myVemap.FindIndex(v2);
1572       if (iv2 == 0) iv2 = myVemap.Add(v2);
1573       iv3 = myVemap.FindIndex(v3);
1574       if (iv3 == 0) iv3 = myVemap.Add(v3);
1575
1576       if (orFace == TopAbs_REVERSED) Tri(i++).Set(iv1, iv3, iv2);
1577       else Tri(i++).Set(iv1, iv2, iv3);
1578     }
1579
1580     Standard_Integer nbVertices = myVemap.Extent();
1581     Handle(Poly_Triangulation) T = new Poly_Triangulation(nbVertices, nTri, Standard_True);
1582     Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
1583     Trian = Tri;
1584     TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
1585     TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
1586
1587     for (i = 1; i <= nbVertices; i++) {
1588       index = myVemap.FindKey(i);
1589       Nodes(i) = Pnt(index);
1590       Nodes2d(i).SetXY(Vertex(index).Coord());
1591     }
1592
1593     T->Deflection(theDefFace);
1594
1595     // storage of triangulation in BRep.
1596     BRep_Builder B1;
1597     //TopLoc_Location loc = theFace.Location();
1598     if (!loc.IsIdentity()) {
1599       gp_Trsf tr = loc.Transformation();
1600       tr.Invert();
1601       for (i = Nodes.Lower(); i <= Nodes.Upper(); i++)
1602         Nodes(i).Transform(tr);
1603     }
1604     B1.UpdateFace(theFace, T);
1605
1606     // implement polygons on triangulation in the face:
1607     BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(myInternaledges);
1608
1609     for (; It.More(); It.Next()) {
1610       const BRepMesh_PairOfPolygon& pair = It.Value();
1611       const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
1612       const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
1613       if ( NOD1 == NOD2 ) {
1614         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly, TOld,loc);
1615         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
1616       }
1617       else {
1618         B.UpdateEdge(TopoDS::Edge(It.Key()), NullPoly,   TOld,loc);
1619         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
1620       }
1621     }
1622   }
1623  }
1624
1625  catch(Standard_Failure)
1626  {
1627    //   MESH_FAILURE(theFace);
1628  }
1629 }
1630
1631
1632 //=======================================================================
1633 //function : Triangle
1634 //purpose  :
1635 //=======================================================================
1636
1637 const BRepMesh_Triangle& BRepMesh_FastDiscretFace::Triangle(const Standard_Integer Index) const
1638 {
1639   return myStructure->GetElement(Index);
1640 }
1641
1642 //=======================================================================
1643 //function : NbEdges
1644 //purpose  :
1645 //=======================================================================
1646
1647 /*Standard_Integer BRepMesh_FastDiscretFace::NbEdges() const
1648 {
1649   return myStructure->NbLinks();
1650 }*/
1651
1652 //=======================================================================
1653 //function : Edge
1654 //purpose  :
1655 //=======================================================================
1656
1657 const BRepMesh_Edge& BRepMesh_FastDiscretFace::Edge(const Standard_Integer Index) const
1658 {
1659   return myStructure->GetLink(Index);
1660 }
1661
1662
1663 //=======================================================================
1664 //function : Vertex
1665 //purpose  :
1666 //=======================================================================
1667
1668 const BRepMesh_Vertex& BRepMesh_FastDiscretFace::Vertex(const Standard_Integer Index) const
1669 {
1670   return myStructure->GetNode(Index);
1671 }
1672
1673 //=======================================================================
1674 //function : Pnt
1675 //purpose  :
1676 //=======================================================================
1677
1678 const gp_Pnt& BRepMesh_FastDiscretFace::Pnt(const Standard_Integer Index) const
1679 {
1680   return myLocation3d(myStructure->GetNode(Index).Location3d());
1681 }
1682
1683 //=======================================================================
1684 //function : FindUV
1685 //purpose  :
1686 //=======================================================================
1687
1688 gp_XY BRepMesh_FastDiscretFace::FindUV(const TopoDS_Vertex&                 theV,
1689                                        const gp_Pnt2d&                      theXY,
1690                                        const Standard_Integer               theIp,
1691                                        const Handle(BRepAdaptor_HSurface)&  theSFace,
1692                                        const Standard_Real                  theMinDist,
1693                                        BRepMesh_DataMapOfIntegerListOfXY&   theLocation2dMap)
1694 {
1695   gp_XY anUV;
1696   if (theLocation2dMap.IsBound(theIp))
1697   {
1698     BRepMesh_ListOfXY& L = theLocation2dMap.ChangeFind(theIp);
1699     anUV = L.First();
1700     if (L.Extent() != 1)
1701     {
1702       BRepMesh_ListIteratorOfListOfXY it(L);
1703       it.Next();
1704       Standard_Real dd, dmin = theXY.Distance(gp_Pnt2d(anUV));
1705       for (; it.More(); it.Next())
1706       {
1707         dd = theXY.Distance(gp_Pnt2d(it.Value()));
1708         if (dd < dmin)
1709         {
1710           anUV = it.Value();
1711           dmin = dd;
1712         }
1713       }
1714     }
1715
1716     const Standard_Real tol = Min(2. * BRep_Tool::Tolerance(theV), theMinDist);
1717
1718     const Standard_Real Utol2d = .5 * (theSFace->LastUParameter() - theSFace->FirstUParameter());
1719     const Standard_Real Vtol2d = .5 * (theSFace->LastVParameter() - theSFace->FirstVParameter());
1720
1721     const gp_Pnt p1 = theSFace->Value(anUV.X(), anUV.Y());
1722     const gp_Pnt p2 = theSFace->Value(theXY.X(), theXY.Y());
1723
1724     if (Abs(anUV.X() - theXY.X()) > Utol2d ||
1725         Abs(anUV.Y() - theXY.Y()) > Vtol2d ||
1726         !p1.IsEqual(p2, tol))
1727     {
1728       anUV = theXY.Coord();
1729       L.Append(anUV);
1730     }
1731   }
1732   else
1733   {
1734     anUV = theXY.Coord();
1735     BRepMesh_ListOfXY L;
1736     L.Append(anUV);
1737     theLocation2dMap.Bind(theIp, L);
1738   }
1739   return anUV;
1740 }
1741
1742
1743 static Standard_Boolean GetVertexParameters(const TopoDS_Vertex& theVert,
1744                                             const TopoDS_Face& theFace,
1745                                             gp_Pnt2d& thePoint)
1746 {
1747   TopLoc_Location L;
1748   const Handle(Geom_Surface)& S = BRep_Tool::Surface(theFace,L);
1749   L = L.Predivided(theVert.Location());
1750   BRep_ListIteratorOfListOfPointRepresentation itpr =
1751     ((*((Handle(BRep_TVertex)*) &theVert.TShape()))->Points());
1752   // Check first if there are PointRepresentation (case non Manifold)
1753
1754   while (itpr.More()) {
1755     if (itpr.Value()->IsPointOnSurface(S,L)) {
1756       thePoint.SetCoord(itpr.Value()->Parameter(),
1757                         itpr.Value()->Parameter2());
1758       return Standard_True;
1759     }
1760     itpr.Next();
1761   }
1762   return Standard_False;
1763 }
1764
1765 //=======================================================================
1766 //function : Add
1767 //purpose  : method intended to add internal myVertices in triangulation.
1768 //=======================================================================
1769 void BRepMesh_FastDiscretFace::Add(const TopoDS_Vertex&                theVert,
1770                                    const TopoDS_Face&                  theFace,
1771                                    const Handle(BRepAdaptor_HSurface)& thegFace)
1772 {
1773   const TopAbs_Orientation anOrient = theVert.Orientation();
1774   gp_Pnt2d uvXY;
1775   if( anOrient != TopAbs_INTERNAL || !GetVertexParameters(theVert,theFace,uvXY))
1776     return;
1777   Standard_Integer indVert =0;
1778   if (myVertices.IsBound(theVert))
1779     indVert = myVertices.Find(theVert);
1780   else
1781   {
1782     myNbLocat++;
1783     myLocation3d.Bind(myNbLocat, BRep_Tool::Pnt(theVert));
1784     indVert = myNbLocat;
1785     myVertices.Bind(theVert, indVert);
1786   }
1787   Standard_Real mindist = BRep_Tool::Tolerance(theVert);
1788   // gp_Pnt2d uvXY = BRep_Tool::Parameters(theVert,theFace);
1789   gp_XY anUV = FindUV(theVert, uvXY, indVert, thegFace, mindist, myLocation2d);
1790   BRepMesh_Vertex vf(anUV, indVert, BRepMesh_Fixed);
1791   Standard_Integer ivff = myStructure->AddNode(vf);
1792   Standard_Integer isvf = myVemap.FindIndex(ivff);
1793   if (isvf == 0) isvf = myVemap.Add(ivff);
1794 }