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