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