src/BRepMesh/BRepMesh_FastDiscretFace.cxx

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