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