src/BRepMesh/BRepMesh_FastDiscret.cxx

   1 // Created on: 1996-02-27
   2 // Created by: Ekaterina SMIRNOVA
   3 // Copyright (c) 1996-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <BRepMesh_FastDiscret.ixx>
  18
  19 #include <BRepMesh_FastDiscretFace.hxx>
  20 #include <BRepMesh_FaceAttribute.hxx>
  21 #include <BRepMesh_DataStructureOfDelaun.hxx>
  22 #include <BRepMesh_ClassifierPtr.hxx>
  23 #include <BRepMesh_Classifier.hxx>
  24 #include <BRepMesh_WireChecker.hxx>
  25 #include <BRepMesh_GeomTool.hxx>
  26 #include <BRepMesh_PairOfPolygon.hxx>
  27 #include <BRepMesh_DataMapOfShapePairOfPolygon.hxx>
  28 #include <BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon.hxx>
  29 #include <Geom_Plane.hxx>
  30 #include <GeomAbs_IsoType.hxx>
  31 #include <GeomAbs_SurfaceType.hxx>
  32 #include <TopAbs.hxx>
  33 #include <TColStd_HArray1OfReal.hxx>
  34 #include <Precision.hxx>
  35
  36 #include <BRep_Builder.hxx>
  37 #include <BRep_Tool.hxx>
  38 #include <Poly_Triangulation.hxx>
  39 #include <Poly_PolygonOnTriangulation.hxx>
  40 #include <Poly_Connect.hxx>
  41 #include <TColStd_SequenceOfInteger.hxx>
  42 #include <TColStd_Array1OfInteger.hxx>
  43 #include <TColStd_HArray1OfInteger.hxx>
  44
  45 #include <TColgp_Array1OfPnt.hxx>
  46 #include <TColgp_Array1OfPnt2d.hxx>
  47 #include <Precision.hxx>
  48
  49 #include <BRepAdaptor_Curve.hxx>
  50 #include <BRepAdaptor_Surface.hxx>
  51 #include <BRepAdaptor_HSurface.hxx>
  52 #include <BRepTools.hxx>
  53 #include <BndLib_Add3dCurve.hxx>
  54 #include <BRepBndLib.hxx>
  55 #include <Bnd_Box.hxx>
  56 #include <TopoDS.hxx>
  57 #include <TopExp.hxx>
  58 #include <TopExp_Explorer.hxx>
  59
  60 #include <Geom2d_Curve.hxx>
  61
  62 #include <TColStd_DataMapOfIntegerInteger.hxx>
  63 #include <BRepMesh_ShapeTool.hxx>
  64 #include <ElSLib.hxx>
  65 #include <Geom_Surface.hxx>
  66 #include <Adaptor3d_IsoCurve.hxx>
  67 #include <BRepMesh_IndexedMapOfVertex.hxx>
  68 #include <Extrema_LocateExtPC.hxx>
  69
  70 #include <BRepMesh_ListOfXY.hxx>
  71 #include <BRepMesh_ListIteratorOfListOfXY.hxx>
  72
  73 #include <TColStd_Array1OfInteger.hxx>
  74 #include <Standard_ErrorHandler.hxx>
  75 #include <Standard_Failure.hxx>
  76 //#include <TColStd_DataMapOfInteger.hxx>
  77 #include <TColGeom2d_SequenceOfCurve.hxx>
  78 #include <TopTools_SequenceOfShape.hxx>
  79 #include <NCollection_IncAllocator.hxx>
  80
  81 #include <BRep_ListIteratorOfListOfPointRepresentation.hxx>
  82 #include <BRep_PointRepresentation.hxx>
  83 #include <BRep_TVertex.hxx>
  84 #include <TColStd_MapOfInteger.hxx>
  85 #include <SortTools_ShellSortOfReal.hxx>
  86 #include <TCollection_CompareOfReal.hxx>
  87
  88 #include <TopTools_HArray1OfShape.hxx>
  89 #include <TopTools_ListIteratorOfListOfShape.hxx>
  90
  91 #include <vector>
  92
  93 #ifdef HAVE_TBB
  94   // paralleling using Intel TBB
  95   #include <tbb/parallel_for_each.h>
  96 #endif
  97
  98 #define UVDEFLECTION 1.e-05
  99
 100 inline Standard_Real MaxFaceTol (const TopoDS_Face& theFace)
 101 {
 102   Standard_Real T, TMax = BRep_Tool::Tolerance(theFace);
 103   TopExp_Explorer Ex;
 104
 105   for (Ex.Init(theFace,TopAbs_EDGE); Ex.More(); Ex.Next())
 106   {
 107     T = BRep_Tool::Tolerance(TopoDS::Edge(Ex.Current()));
 108     if (T > TMax) TMax = T;
 109   }
 110
 111   for (Ex.Init(theFace,TopAbs_VERTEX); Ex.More(); Ex.Next())
 112   {
 113     T = BRep_Tool::Tolerance(TopoDS::Vertex(Ex.Current()));
 114     if (T > TMax) TMax = T;
 115   }
 116
 117   return TMax;
 118 }
 119
 120
 121 //=======================================================================
 122 //function : BRepMesh_FastDiscret
 123 //purpose  :
 124 //=======================================================================
 125 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    theDefle,
 126                                            const Standard_Real    theAngl,
 127                                            const Bnd_Box&         theBox,
 128                                            const Standard_Boolean theWithShare,
 129                                            const Standard_Boolean theInshape,
 130                                            const Standard_Boolean theRelative,
 131                                            const Standard_Boolean theShapetrigu,
 132                                            const Standard_Boolean isInParallel)
 133 : myAngle (theAngl),
 134   myDeflection (theDefle),
 135   myWithShare (theWithShare),
 136   myInParallel (isInParallel),
 137   myNbLocat (0),
 138   myRelative (theRelative),
 139   myShapetrigu (theShapetrigu),
 140   myInshape (theInshape)
 141 {
 142   myAllocator = new NCollection_IncAllocator(64000);
 143   if(myRelative)
 144     BoxMaxDimension(theBox, myDtotale);
 145 }
 146
 147 //=======================================================================
 148 //function : BRepMesh_FastDiscret
 149 //purpose  :
 150 //=======================================================================
 151
 152 BRepMesh_FastDiscret::BRepMesh_FastDiscret(const Standard_Real    theDefle,
 153                                            const TopoDS_Shape&    theShape,
 154                                            const Bnd_Box&         theBox,
 155                                            const Standard_Real    theAngl,
 156                                            const Standard_Boolean theWithShare,
 157                                            const Standard_Boolean theInshape,
 158                                            const Standard_Boolean theRelative,
 159                                            const Standard_Boolean theShapetrigu,
 160                                            const Standard_Boolean isInParallel)
 161 : myAngle (theAngl),
 162   myDeflection (theDefle),
 163   myWithShare (theWithShare),
 164   myInParallel (isInParallel),
 165   myNbLocat (0),
 166   myRelative (theRelative),
 167   myShapetrigu (theShapetrigu),
 168   myInshape (theInshape)
 169 {
 170   myAllocator = new NCollection_IncAllocator(64000);
 171   if(myRelative)
 172     BoxMaxDimension(theBox, myDtotale);
 173   Perform(theShape);
 174 }
 175
 176 //=======================================================================
 177 //function : SetParallel
 178 //purpose  :
 179 //=======================================================================
 180 void BRepMesh_FastDiscret::SetParallel (const Standard_Boolean theInParallel)
 181 {
 182   myInParallel = theInParallel;
 183 }
 184
 185 //=======================================================================
 186 //function : IsParallel
 187 //purpose  :
 188 //=======================================================================
 189 Standard_Boolean BRepMesh_FastDiscret::IsParallel() const
 190 {
 191   return myInParallel;
 192 }
 193
 194 //=======================================================================
 195 //function : BoxMaxDimension
 196 //purpose  :
 197 //=======================================================================
 198
 199 void BRepMesh_FastDiscret::BoxMaxDimension(const Bnd_Box& theBox, Standard_Real& theMaxDim)
 200 {
 201   if(theBox.IsVoid())
 202     return;
 203   Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
 204   theBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
 205   theMaxDim = TXmax-TXmin;
 206   const Standard_Real dy = TYmax-TYmin;
 207   const Standard_Real dz = TZmax-TZmin;
 208   if (dy > theMaxDim) theMaxDim = dy;
 209   if (dz > theMaxDim) theMaxDim = dz;
 210 }
 211
 212 //=======================================================================
 213 //function : RelativeEdgeDeflection
 214 //purpose  :
 215 //=======================================================================
 216
 217 Standard_Real BRepMesh_FastDiscret::RelativeEdgeDeflection(const TopoDS_Edge& theEdge,
 218                                                                                                    const Standard_Real theDefle,
 219                                                                                                      const Standard_Real theDTotale,
 220                                                            Standard_Real& theDefCoef)
 221 {
 222   theDefCoef = 1.;
 223   Standard_Real defedge = theDefle;
 224   if(theEdge.IsNull())
 225     return defedge;
 226
 227   Bnd_Box B;
 228   BRepBndLib::Add(theEdge, B);
 229   BoxMaxDimension(B, defedge);
 230
 231   // adjusted in relation to the total size:
 232   theDefCoef = theDTotale/(2*defedge);
 233   if (theDefCoef < 0.5) theDefCoef = 0.5;
 234   if (theDefCoef > 2.) theDefCoef = 2.;
 235   defedge = theDefCoef * defedge * theDefle;
 236
 237   return defedge;
 238 }
 239
 240 //=======================================================================
 241 //function : Perform(shape)
 242 //purpose  :
 243 //=======================================================================
 244
 245 void BRepMesh_FastDiscret::Perform(const TopoDS_Shape& theShape)
 246 {
 247   TopTools_IndexedDataMapOfShapeListOfShape anAncestors;
 248   TopExp::MapShapesAndAncestors(theShape, TopAbs_EDGE, TopAbs_FACE, anAncestors);
 249   std::vector<TopoDS_Face> aFaces;
 250   for (TopExp_Explorer ex(theShape, TopAbs_FACE); ex.More(); ex.Next()) {
 251     TopoDS_Face aF = TopoDS::Face(ex.Current());
 252     Add(aF, anAncestors);
 253     aFaces.push_back(aF);
 254   }
 255
 256   if (myInParallel)
 257   {
 258   #ifdef HAVE_TBB
 259     CreateMutexesForSubShapes(theShape, TopAbs_EDGE);
 260     // mesh faces in parallel threads using TBB
 261     tbb::parallel_for_each (aFaces.begin(), aFaces.end(), *this);
 262   #else
 263     // alternative parallelization not yet available
 264     for (std::vector<TopoDS_Face>::iterator it(aFaces.begin()); it != aFaces.end(); it++)
 265       Process (*it);
 266   #endif
 267     RemoveAllMutexes();
 268   }
 269   else
 270   {
 271     for (std::vector<TopoDS_Face>::iterator it(aFaces.begin()); it != aFaces.end(); it++)
 272       Process (*it);
 273   }
 274 }
 275
 276
 277 //=======================================================================
 278 //function : Process
 279 //purpose  :
 280 //=======================================================================
 281
 282 void BRepMesh_FastDiscret::Process(const TopoDS_Face& theFace) const
 283 {
 284   //cout << "START face " << theFace.TShape().operator->() << endl << flush;
 285   Handle(BRepMesh_FaceAttribute) fattribute;
 286   if ( GetFaceAttribute (theFace, fattribute) )
 287   {
 288     BRepMesh_FastDiscretFace aTool (GetAngle(), WithShare());
 289     aTool.Add (theFace, fattribute, GetMapOfDefEdge(), myMutexProvider);
 290   }
 291   //cout << "END   face " << theFace.TShape().operator->() << endl << flush;
 292 }
 293
 294 //=======================================================================
 295 //function : Add(face)
 296 //purpose  :
 297 //=======================================================================
 298
 299 #define MESH_FAILURE(theface)        \
 300 {                                    \
 301   myFacestate = BRepMesh_Failure;    \
 302   myNottriangulated.Append(theface); \
 303   return;                            \
 304 }
 305
 306 void BRepMesh_FastDiscret::Add(const TopoDS_Face& theface,
 307                                const TopTools_IndexedDataMapOfShapeListOfShape& theAncestors)
 308 {
 309 #ifndef DEB_MESH
 310   try
 311   {
 312     OCC_CATCH_SIGNALS
 313 #endif
 314   TopoDS_Face face = theface;
 315   BRepTools::Update(face);
 316   face.Orientation(TopAbs_FORWARD);
 317   myStructure.Nullify();
 318   Handle(NCollection_IncAllocator) anAlloc = Handle(NCollection_IncAllocator)::DownCast(myAllocator);
 319   anAlloc->Reset(Standard_False);
 320   myStructure=new BRepMesh_DataStructureOfDelaun(anAlloc);
 321
 322   Standard_Real aUmin, aVmin, aUmax, aVmax;
 323   BRepTools::UVBounds (theface, aUmin, aUmax, aVmin, aVmax);
 324   Standard_Real aTolU = Max( Precision::PConfusion(), (aUmax - aUmin) * UVDEFLECTION );
 325   Standard_Real aTolV = Max( Precision::PConfusion(), (aVmax - aVmin) * UVDEFLECTION );
 326   myStructure->Data().SetCellSize ( 14 * aTolU, 14 * aTolV );
 327   myStructure->Data().SetTolerance( aTolU, aTolV );
 328
 329   BRepAdaptor_Surface  BS(face, Standard_False);
 330   Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
 331
 332   GeomAbs_SurfaceType thetype;
 333   thetype = BS.GetType();
 334
 335   gp_Pnt2d uvFirst, uvLast;
 336
 337   Handle(Poly_Triangulation) T;
 338   TopLoc_Location loc;
 339
 340   if (!myWithShare) {
 341     myVertices.Clear();
 342     myEdges.Clear();
 343   }
 344
 345   myVemap.Clear();
 346   myLocation2d.Clear();
 347   myInternaledges.Clear();
 348
 349   Standard_Integer i;
 350   i = 1;
 351
 352   Standard_Real defedge, defface;
 353   Standard_Integer nbEdge = 0;
 354   Standard_Real savangle = myAngle;
 355   Standard_Real cdef;
 356   Standard_Real maxdef = 2.* MaxFaceTol(theface);
 357   defface = 0.;
 358
 359   if (!myRelative) defface = Max(myDeflection, maxdef);
 360
 361   TColStd_SequenceOfReal aFSeq, aLSeq;
 362   TColGeom2d_SequenceOfCurve aCSeq;
 363   TopTools_SequenceOfShape aShSeq;
 364
 365   TopoDS_Iterator exW(face);
 366
 367   for (; exW.More(); exW.Next()) {
 368     const TopoDS_Shape& aWire = exW.Value();
 369     if (aWire.ShapeType() != TopAbs_WIRE)
 370       continue;
 371     TopoDS_Iterator ex(aWire);
 372     for(; ex.More(); ex.Next()) {
 373       const TopoDS_Edge& edge = TopoDS::Edge(ex.Value());
 374       nbEdge++;
 375       if (!myMapdefle.IsBound(edge)) {
 376         if (myRelative) {
 377           if (myEdges.IsBound(edge)) {
 378             const BRepMesh_PairOfPolygon& pair = myEdges.Find(edge);
 379             const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
 380             defedge = P->Deflection();
 381           }
 382           else {
 383             defedge = RelativeEdgeDeflection(edge, myDeflection, myDtotale, cdef);
 384             myAngle = savangle * cdef;
 385           }
 386           defface = defface + defedge;
 387           defface = Max(maxdef, defface);
 388         }
 389         else defedge = myDeflection;
 390
 391         defedge = Max(maxdef, defedge);
 392         defedge = Max(UVDEFLECTION , defedge);
 393         myMapdefle.Bind(edge, defedge);
 394       }
 395       else{
 396         defedge = myMapdefle(edge);
 397         if (myRelative) {defface = defface + defedge; defface = Max(maxdef, defface);}
 398       }
 399       Standard_Real f1,l1;
 400       Handle(Geom2d_Curve) C = BRep_Tool::CurveOnSurface(edge, face, f1, l1);
 401       if (C.IsNull()) continue;
 402
 403       aFSeq.Append(f1);
 404       aLSeq.Append(l1);
 405       aCSeq.Append(C);
 406       aShSeq.Append(edge);
 407       Add(edge, face, gFace, C, theAncestors, defedge, f1, l1);
 408       myAngle = savangle;
 409     }
 410   }
 411
 412   if (nbEdge == 0 || myVemap.Extent() < 3)
 413   {
 414       MESH_FAILURE(theface);
 415   }
 416
 417   if (myRelative ) defface = defface / nbEdge;
 418   else             defface = myDeflection;
 419
 420   if (myWithShare) defface = Max(maxdef, defface);
 421
 422   T = BRep_Tool::Triangulation(face, loc);
 423
 424   if (!myShapetrigu || T.IsNull()) {
 425
 426     Standard_Real xCur, yCur;
 427     Standard_Real maxX, minX, maxY, minY;
 428     minX=minY=1.e100;
 429     maxX=maxY=-1.e100;
 430
 431     Standard_Integer ipn = 0;
 432     Standard_Integer i1 =1;
 433     for (i1 = 1; i1 <= myVemap.Extent(); i1++) {
 434       const BRepMesh_Vertex& aV = myStructure->GetNode(myVemap.FindKey(i1));
 435       ipn++;
 436       xCur=aV.Coord().X();
 437       yCur=aV.Coord().Y();
 438       minX=Min(xCur, minX);
 439       maxX=Max(xCur, maxX);
 440       minY=Min(yCur, minY);
 441       maxY=Max(yCur, maxY);
 442     }
 443     Standard_Real myumin = minX;
 444     Standard_Real myumax = maxX;
 445     Standard_Real myvmin = minY;
 446     Standard_Real myvmax = maxY;
 447
 448     const Standard_Real umin = BS.FirstUParameter();
 449     const Standard_Real umax = BS.LastUParameter();
 450     const Standard_Real vmin = BS.FirstVParameter();
 451     const Standard_Real vmax = BS.LastVParameter();
 452
 453     if (myumin < umin || myumax > umax)
 454     {
 455       if (BS.IsUPeriodic())
 456       {
 457         if ((myumax - myumin) > (umax - umin))
 458         {
 459           myumax = myumin + (umax - umin);
 460         }
 461       }
 462       else
 463       {
 464         if (umin > myumin) myumin = umin;
 465         if (umax < myumax) myumax = umax;
 466       }
 467     }
 468
 469     if (myvmin < vmin || myvmax > vmax)
 470     {
 471       if (BS.IsVPeriodic())
 472       {
 473         if ((myvmax - myvmin) > (vmax - vmin))
 474         {
 475           myvmax = myvmin + (vmax - vmin);
 476         }
 477       }
 478       else
 479       {
 480         if (vmin > myvmin) myvmin = vmin;
 481         if (vmax < myvmax) myvmax = vmax;
 482       }
 483     }
 484
 485     // fast verification of the validity of calculated limits. If wrong,
 486     // sure a problem of pcurve.
 487     if (thetype == GeomAbs_BezierSurface &&
 488         (myumin < -0.5 || myumax > 1.5 || myvmin < -0.5 || myvmax > 1.5))
 489     {
 490       MESH_FAILURE(theface);
 491     }
 492
 493     //define parameters for correct parametrics
 494
 495     Standard_Real deltaX = 1.0;
 496     Standard_Real deltaY = 1.0;
 497     Standard_Integer nbVertices = myVemap.Extent();
 498     const Standard_Real tolclass = Precision::PConfusion(); //0.03*Max(myumax-myumin, myvmax-myvmin);
 499
 500     BRepMesh_ClassifierPtr classifier = new BRepMesh_Classifier;
 501     {
 502       BRepMesh_WireChecker aDFaceChecker(face,
 503         tolclass, myInternaledges, myVemap, myStructure,
 504         myumin, myumax, myvmin, myvmax, myInParallel);
 505       aDFaceChecker.ReCompute(classifier);
 506
 507       myFacestate = aDFaceChecker.Status();
 508       if (myFacestate == BRepMesh_SelfIntersectingWire)
 509       {
 510         Standard_Integer nbmaill = 0;
 511         Standard_Real eps = Precision::Confusion();
 512         while (nbmaill < 5 && myFacestate != BRepMesh_ReMesh)
 513         {
 514           nbmaill++;
 515
 516           //clear the structure of links
 517           myStructure.Nullify();
 518           myStructure = new BRepMesh_DataStructureOfDelaun(anAlloc);
 519
 520           myVemap.Clear();
 521           myLocation2d.Clear();
 522           myInternaledges.Clear();
 523
 524           Standard_Integer j1;
 525           for(j1 = 1; j1 <= aShSeq.Length(); j1++)
 526           {
 527             const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
 528             if (myEdges.IsBound(edge))
 529             {
 530               myEdges.UnBind(edge);
 531               myInternaledges.UnBind(edge);
 532             }
 533           }
 534
 535
 536           for( j1 = 1; j1 <= aShSeq.Length(); j1++)
 537           {
 538             const TopoDS_Edge& edge = TopoDS::Edge(aShSeq.Value(j1));
 539             defedge = myMapdefle(edge) / 3.;
 540             defedge = Max(defedge, eps);
 541             myMapdefle.Bind(edge, defedge);
 542             const Handle(Geom2d_Curve)& C = aCSeq.Value(j1);
 543             Add(edge, face, gFace, C, theAncestors, defedge, aFSeq.Value(j1), aLSeq.Value(j1));
 544           }
 545
 546           aDFaceChecker.ReCompute(classifier);
 547           if (aDFaceChecker.Status() == BRepMesh_NoError)
 548           {
 549             myFacestate = BRepMesh_ReMesh;
 550           }
 551           nbVertices = myVemap.Extent();
 552         }
 553       }
 554     }
 555
 556     if (myFacestate != BRepMesh_NoError && myFacestate != BRepMesh_ReMesh)
 557     {
 558       myNottriangulated.Append(face);
 559       classifier.Nullify();
 560       return;
 561     }
 562
 563     // try to find the real length:
 564     // akm (bug OCC16) : We must calculate these measures in non-singular
 565     //     parts of face. Let's try to compute average value of three
 566     //     (umin, (umin+umax)/2, umax), and respectively for v.
 567     //                 vvvvv
 568     Standard_Real longu = 0.0, longv = 0.0; //, last , first;
 569     gp_Pnt P11, P12, P21, P22, P31, P32;
 570
 571     Standard_Real du = (myumax-myumin)/20;
 572     Standard_Real dv = (myvmax-myvmin)/20;
 573     Standard_Real dfuave=(myumin+myumax)/2, dfucur;
 574     Standard_Real dfvave=(myvmin+myvmax)/2, dfvcur;
 575     // U loop
 576     BS.D0 (myumin, myvmin, P11);
 577     BS.D0 (myumin, dfvave, P21);
 578     BS.D0 (myumin, myvmax, P31);
 579     for (i1=1, dfucur=myumin+du; i1 <= 20; i1++, dfucur+=du)  {
 580       BS.D0 (dfucur, myvmin, P12);
 581       BS.D0 (dfucur, dfvave, P22);
 582       BS.D0 (dfucur, myvmax, P32);
 583       longu += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
 584       P11 = P12;
 585       P21 = P22;
 586       P31 = P32;
 587     }
 588     // V loop
 589     BS.D0(myumin, myvmin, P11);
 590     BS.D0(dfuave, myvmin, P21);
 591     BS.D0(myumax, myvmin, P31);
 592     for (i1=1, dfvcur=myvmin+dv; i1 <= 20; i1++, dfvcur+=dv) {
 593       BS.D0 (myumin, dfvcur, P12);
 594       BS.D0 (dfuave, dfvcur, P22);
 595       BS.D0 (myumax, dfvcur, P32);
 596       longv += ( P11.Distance(P12) + P21.Distance(P22) + P31.Distance(P32) );
 597       P11 = P12;
 598       P21 = P22;
 599       P31 = P32;
 600     }
 601     longu /= 3.;
 602     longv /= 3.;
 603     // akm (bug OCC16) ^^^^^
 604
 605     if (longu <= 1.e-16 || longv <= 1.e-16) {
 606       //yes, it is seen!!
 607 #ifdef DEB_MESH_CHRONO
 608       chMaillEdges.Stop();
 609       MESH_CHRONO_TSTOP(thetype);
 610 #endif
 611       MESH_FAILURE(theface);
 612     }
 613
 614
 615     if (thetype == GeomAbs_Torus)  {
 616       gp_Torus Tor = BS.Torus();
 617       Standard_Real r = Tor.MinorRadius(), R = Tor.MajorRadius();
 618       Standard_Real Du, Dv;//, pasu, pasv;
 619
 620       Dv = Max(1.0e0 - (defface/r),0.0e0) ;
 621       Standard_Real oldDv = 2.0 * ACos (Dv);
 622       oldDv = Min(oldDv, myAngle);
 623       Dv  =  0.9*oldDv; //TWOTHIRD * oldDv;
 624       Dv = oldDv;
 625
 626       Standard_Integer nbV = Max((Standard_Integer)((myvmax-myvmin)/Dv), 2);
 627       Dv = (myvmax-myvmin)/(nbV+1);
 628
 629       Standard_Real ru = R + r;
 630       if (ru > 1.e-16) {
 631         Du = Max(1.0e0 - (defface/ru),0.0e0);
 632         Du  = (2.0 * ACos (Du));
 633         Du = Min(Du, myAngle);
 634         Standard_Real aa = sqrt(Du*Du + oldDv*oldDv);
 635         if(aa < gp::Resolution())
 636           return;
 637
 638         Du = Du * Min(oldDv, Du) / aa;
 639       }
 640       else Du = Dv;
 641
 642       Standard_Integer nbU = Max((Standard_Integer)((myumax-myumin)/Du), 2);
 643       nbU = Max(nbU, (Standard_Integer)(nbV*(myumax-myumin)*R/((myvmax-myvmin)*r)/5.));
 644
 645       Du = (myumax-myumin)/(nbU+1);
 646       //-- DeltaX and DeltaY are chosen so that to avoid "jumping"
 647       //-- of points on the grid
 648       deltaX = Du;
 649       deltaY = Dv;
 650     }
 651     else if (thetype == GeomAbs_Cylinder) {
 652       /*Standard_Real aMax = Max(myumax,myvmax);
 653       deltaX = 0.01;  //1./aMax; //0.01;
 654       deltaY = 1.0;*/
 655       gp_Cylinder Cyl = BS.Cylinder();
 656       Standard_Real R = Cyl.Radius();
 657       // Calculate parameters for iteration in U direction
 658       Standard_Real Du = 1.0 - (defface/R);
 659       if (Du < 0.0) Du = 0.0;
 660       Du = 2.0 * ACos (Du);
 661       if (Du > GetAngle()) Du = GetAngle();
 662       deltaX = Du/longv;
 663       deltaY = 1.;
 664     }
 665     else {
 666       deltaX = (myumax-myumin)/longu;
 667       deltaY = (myvmax-myvmin)/longv;
 668     }
 669
 670     if(!myMapattrib.IsBound(theface))
 671     {
 672       Handle(BRepMesh_FaceAttribute) aFA = new BRepMesh_FaceAttribute();
 673       aFA->GetDefFace() = defface;
 674       aFA->GetUMax() = myumax;
 675       aFA->GetVMax() = myvmax;
 676       aFA->GetUMin() = myumin;
 677       aFA->GetVMin() = myvmin;
 678       aFA->GetDeltaX() = deltaX;
 679       aFA->GetDeltaY() = deltaY;
 680       aFA->GetMinX() = minX;
 681       aFA->GetMinY() = minY;
 682       aFA->GetClassifier() = classifier;
 683       myMapattrib.Bind(theface, aFA);
 684     }
 685
 686     //Standard_Integer nbVertices = myVemap.Extent();
 687     T = new Poly_Triangulation(nbVertices, 1, Standard_True);
 688     TColgp_Array1OfPnt&  Nodes = T->ChangeNodes();
 689     TColgp_Array1OfPnt2d& Nodes2d = T->ChangeUVNodes();
 690
 691     Standard_Integer index;
 692     for (i = 1; i <= nbVertices; i++) {
 693       index = myVemap.FindKey(i);
 694       Nodes(i) = Pnt(index);
 695       Nodes2d(i).SetXY(Vertex(index).Coord());
 696     }
 697
 698     // storage of triangulation in the BRep.
 699     //TopLoc_Location loc = face.Location();
 700     if (!loc.IsIdentity()) {
 701       gp_Trsf tr = loc.Transformation();
 702       tr.Invert();
 703       for (i = Nodes.Lower(); i <= Nodes.Upper(); i++)
 704         Nodes(i).Transform(tr);
 705     }
 706
 707     BRep_Builder B;
 708     B.UpdateFace(face, T);
 709
 710     BRepMesh_DataMapIteratorOfDataMapOfShapePairOfPolygon It(myInternaledges);
 711     for (; It.More(); It.Next()) {
 712       const BRepMesh_PairOfPolygon& pair = It.Value();
 713       const Handle(Poly_PolygonOnTriangulation)& NOD1 = pair.First();
 714       const Handle(Poly_PolygonOnTriangulation)& NOD2 = pair.Last();
 715       if ( NOD1 == NOD2 )
 716         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, T, loc);
 717       else
 718         B.UpdateEdge(TopoDS::Edge(It.Key()), NOD1, NOD2, T, loc);
 719     }
 720   }
 721
 722 #ifndef DEB_MESH
 723   }
 724   catch(Standard_Failure)
 725   {
 726     MESH_FAILURE(theface);
 727   }
 728 #endif // DEB_MESH
 729   myStructure.Nullify();
 730 }
 731
 732 //=======================================================================
 733 //function : splitSegment
 734 //purpose  :
 735 //=======================================================================
 736 static void splitSegment( BRepMesh_GeomTool&    theGT,
 737                           const Handle(Geom_Surface)& theSurf,
 738                           const Handle(Geom2d_Curve)& theCurve2d,
 739                           const BRepAdaptor_Curve&    theBAC,
 740                           const Standard_Real         theSquareEDef,
 741                           const Standard_Real         theFirst,
 742                           const Standard_Real         theLast,
 743                           const Standard_Integer      theNbIter)
 744 {
 745   //limit ineration depth
 746   if(theNbIter > 10)
 747     return;
 748   gp_Pnt2d uvf, uvl, uvm;
 749   gp_Pnt   P3dF, P3dL, midP3d, midP3dFromSurf;
 750   Standard_Real midpar;
 751
 752   if(Abs(theLast - theFirst) < 2*Precision::PConfusion())
 753     return;
 754
 755   theCurve2d->D0(theFirst, uvf);
 756   theCurve2d->D0(theLast, uvl);
 757
 758   P3dF = theSurf->Value(uvf.X(), uvf.Y());
 759   P3dL = theSurf->Value(uvl.X(), uvl.Y());
 760
 761   if(P3dF.SquareDistance(P3dL) < theSquareEDef)
 762     return;
 763
 764   uvm = gp_Pnt2d((uvf.XY() + uvl.XY())*0.5);
 765   midP3dFromSurf = theSurf->Value(uvm.X(), uvm.Y());
 766
 767   gp_XYZ aVec = P3dL.XYZ()-P3dF.XYZ();
 768   aVec.Normalize();
 769
 770   gp_XYZ Vec1 = midP3dFromSurf.XYZ() - P3dF.XYZ();
 771   Standard_Real aModulus = Vec1.Dot(aVec);
 772   gp_XYZ aProj = aVec*aModulus;
 773   gp_XYZ aDist = Vec1 - aProj;
 774
 775   if(aDist.SquareModulus() < theSquareEDef)
 776     return;
 777
 778   midpar = (theFirst + theLast) * 0.5;
 779   theBAC.D0(midpar, midP3d);
 780   theGT.AddPoint(midP3d, midpar, Standard_False);
 781
 782   splitSegment(theGT, theSurf, theCurve2d, theBAC, theSquareEDef, theFirst, midpar, theNbIter+1);
 783   splitSegment(theGT, theSurf, theCurve2d, theBAC, theSquareEDef, midpar, theLast,  theNbIter+1);
 784 }
 785
 786 //=======================================================================
 787 //function : Add
 788 //purpose  :
 789 //=======================================================================
 790 void BRepMesh_FastDiscret::Add( const TopoDS_Edge&                  theEdge,
 791                                 const TopoDS_Face&                  theFace,
 792                                 const Handle(BRepAdaptor_HSurface)& theGFace,
 793                                 const Handle(Geom2d_Curve)&         theC2d,
 794                                 const TopTools_IndexedDataMapOfShapeListOfShape& theAncestors,
 795                                 const Standard_Real                 theDefEdge,
 796                                 const Standard_Real                 theFirst,
 797                                 const Standard_Real                 theLast)
 798 {
 799   const TopAbs_Orientation orEdge = theEdge.Orientation();
 800   if (orEdge == TopAbs_EXTERNAL) return;
 801
 802   const Standard_Boolean isEdgeBound = myEdges.IsBound(theEdge);
 803
 804   if (!isEdgeBound)
 805   {
 806     if (Update(theEdge, theFace, theC2d, theDefEdge, theFirst, theLast))
 807     {
 808       return;
 809     }
 810   }
 811
 812   TopoDS_Vertex pBegin, pEnd;
 813   TopExp::Vertices(theEdge, pBegin, pEnd);
 814   if (pBegin.IsNull() || pEnd.IsNull())
 815   {
 816     return;
 817   }
 818
 819   Standard_Real wFirst, wLast;
 820   BRep_Tool::Range(theEdge, theFace, wFirst, wLast);
 821
 822   gp_Pnt2d uvFirst, uvLast;
 823   BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
 824
 825   const Standard_Boolean sameUV = uvFirst.IsEqual(uvLast, Precision::PConfusion());
 826
 827   //Control vertexes tolerances
 828   gp_Pnt pFirst = theGFace->Value(uvFirst.X(), uvFirst.Y());
 829   gp_Pnt pLast  = theGFace->Value(uvLast.X(), uvLast.Y());
 830
 831   Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)),
 832                                     pLast.Distance(BRep_Tool::Pnt(pEnd)));
 833
 834   if(mindist < BRep_Tool::Tolerance(pBegin) ||
 835      mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
 836
 837   // control of degenerated non-coded edges.
 838
 839   Standard_Boolean degener = BRep_Tool::Degenerated(theEdge);
 840
 841   if (!degener)
 842   {
 843     if (pBegin.IsSame(pEnd))
 844     {
 845       // calculation of the length of the edge in 3D
 846       Standard_Real longueur = 0.0;
 847       Standard_Real du = (wLast-wFirst)/20;
 848       gp_Pnt P1, P2;
 849       BRepAdaptor_Curve BC(theEdge);
 850       BC.D0(wFirst, P1);
 851       Standard_Real tolV = BRep_Tool::Tolerance(pBegin);
 852       Standard_Real tolV2 = 1.2*tolV;
 853       for (Standard_Integer l = 1; l <= 20; l++) {
 854         BC.D0(wFirst + l*du, P2);
 855         longueur += P1.Distance(P2);
 856         if (longueur > tolV2) break;
 857         P1 = P2;
 858       }
 859
 860       if (longueur < tolV2)
 861       {
 862         degener = Standard_True;
 863       }
 864     }
 865   }
 866
 867   // Correct UV points
 868   if (sameUV)
 869   {
 870     // 1. is it really sameUV without being degenerated
 871     gp_Pnt2d uvF, uvL;
 872     theC2d->D0(theFirst, uvF);
 873     theC2d->D0(theLast, uvL);
 874     if (!uvFirst.IsEqual(uvF, Precision::PConfusion()))
 875     {
 876       uvFirst = uvF;
 877     }
 878     if (!uvLast.IsEqual(uvL, Precision::PConfusion()))
 879     {
 880       uvLast = uvL;
 881     }
 882   }
 883
 884   gp_XY theUV;
 885
 886   // Process first vertex
 887   Standard_Integer ipf;
 888   if (myVertices.IsBound(pBegin))
 889   {
 890     ipf = myVertices.Find(pBegin);
 891   }
 892   else
 893   {
 894     if (sameUV && myVertices.IsBound(pEnd))
 895     {
 896       ipf = myVertices.Find(pEnd);
 897     }
 898     else
 899     {
 900       myNbLocat++;
 901       ipf = myNbLocat;
 902       myLocation3d.Bind(ipf, BRep_Tool::Pnt(pBegin));
 903     }
 904     myVertices.Bind(pBegin, ipf);
 905   }
 906
 907   Handle(BRepMesh_FaceAttribute) aFaceAttribute;
 908   GetFaceAttribute ( theFace, aFaceAttribute );
 909   theUV = BRepMesh_FastDiscretFace::FindUV(pBegin, uvFirst,
 910     ipf, theGFace, mindist, aFaceAttribute, myLocation2d);
 911
 912   BRepMesh_Vertex vf(theUV, ipf, BRepMesh_Frontier);
 913   Standard_Integer ivf = myStructure->AddNode(vf);
 914
 915   // Process last vertex
 916   Standard_Integer ipl;
 917   if (pEnd.IsSame(pBegin))
 918   {
 919     ipl = ipf;
 920   }
 921   else
 922   {
 923     if (myVertices.IsBound(pEnd))
 924     {
 925       ipl = myVertices.Find(pEnd);
 926     }
 927     else
 928     {
 929       if (sameUV)
 930       {
 931         ipl = ipf;
 932       }
 933       else
 934       {
 935         myNbLocat++;
 936         ipl = myNbLocat;
 937         myLocation3d.Bind(ipl, BRep_Tool::Pnt(pEnd));
 938       }
 939       myVertices.Bind(pEnd,ipl);
 940     }
 941   }
 942
 943   theUV = BRepMesh_FastDiscretFace::FindUV(pEnd, uvLast, ipl,
 944     theGFace, mindist, aFaceAttribute, myLocation2d);
 945
 946   BRepMesh_Vertex vl(theUV, ipl, BRepMesh_Frontier);
 947   Standard_Integer ivl= myStructure->AddNode(vl);
 948
 949   Standard_Integer isvf = myVemap.FindIndex(ivf);
 950   if (isvf == 0) isvf = myVemap.Add(ivf);
 951   Standard_Integer isvl = myVemap.FindIndex(ivl);
 952   if (isvl == 0) isvl = myVemap.Add(ivl);
 953
 954   Standard_Real otherdefedge = 0.5*theDefEdge;
 955   gp_Pnt2d uv;
 956   BRepMesh_Vertex v2;
 957   gp_Pnt   P3d;
 958
 959   Handle(Poly_PolygonOnTriangulation) P1;
 960
 961   if (!isEdgeBound)
 962   {
 963     Handle(Poly_PolygonOnTriangulation) P2;
 964
 965     if (degener)
 966     {
 967       // creation anew:
 968       TColStd_Array1OfInteger Nodes(1, 2), NodInStruct(1, 2);
 969       TColStd_Array1OfReal Param(1, 2);
 970
 971       NodInStruct(1) = ipf;
 972       Nodes(1) = isvf;
 973       Param(1) = wFirst;
 974
 975       NodInStruct(2) = ipl;
 976       Nodes(2) = isvl;
 977       Param(2) = wLast;
 978
 979       P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
 980       P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
 981     }
 982     else
 983     {
 984       if (orEdge == TopAbs_INTERNAL) otherdefedge *= 0.5;
 985
 986       BRepAdaptor_Curve cons;
 987       Standard_Boolean isSameParam = BRep_Tool::SameParameter(theEdge);
 988       if (isSameParam)
 989       {
 990         cons.Initialize(theEdge);
 991       }
 992       else
 993       {
 994         cons.Initialize(theEdge, theFace);
 995       }
 996
 997       TopLoc_Location L;
 998       Standard_Integer nbpmin = 2;
 999       const GeomAbs_CurveType aCurveType = cons.GetType();
1000       if ( aCurveType == GeomAbs_Circle )
1001         nbpmin = 4; //OCC287
1002
1003       BRepMesh_GeomTool GT(cons, wFirst, wLast, 0.5 * myAngle, otherdefedge, nbpmin);
1004
1005       if ( aCurveType == GeomAbs_BSplineCurve )
1006       {
1007         const Standard_Integer aNbInt = cons.NbIntervals( GeomAbs_C1 );
1008         if ( aNbInt > 0 )
1009         {
1010           TColStd_Array1OfReal anIntervals( 1, aNbInt + 1 );
1011           cons.Intervals( anIntervals, GeomAbs_C1 );
1012           for ( Standard_Integer aIntIt = 1; aIntIt <= aNbInt; ++aIntIt )
1013           {
1014             const Standard_Real& aStartInt = anIntervals.Value( aIntIt );
1015             const Standard_Real& anEndInt  = anIntervals.Value( aIntIt + 1 );
1016
1017             BRepMesh_GeomTool aDetalizator( cons, aStartInt, anEndInt,
1018               0.5 * myAngle, otherdefedge, nbpmin );
1019
1020             Standard_Integer aNbAddNodes = aDetalizator.NbPoints();
1021             for ( Standard_Integer aNodeIt = 1; aNodeIt <= aNbAddNodes; ++aNodeIt )
1022             {
1023               Standard_Real aParam;
1024               gp_Pnt        aPoint3d;
1025               gp_Pnt2d      aPoint2d;
1026               aDetalizator.Value( cons, theGFace, aNodeIt, aParam, aPoint3d, aPoint2d );
1027               GT.AddPoint( aPoint3d, aParam, Standard_False );
1028             }
1029           }
1030         }
1031       }
1032
1033       // PTv, chl/922/G9, Take into account internal vertices
1034       // it is necessary for internal edges, which do not split other edges, by their vertex
1035       TopoDS_Iterator exV(theEdge);
1036       for ( ; exV.More(); exV.Next() )
1037       {
1038         TopoDS_Vertex aIntV = TopoDS::Vertex(exV.Value());
1039         if ( aIntV.Orientation() == TopAbs_INTERNAL )
1040          GT.AddPoint(BRep_Tool::Pnt(aIntV),
1041                      BRep_Tool::Parameter(aIntV, theEdge),
1042                      Standard_True);
1043       }
1044
1045       Standard_Integer i;
1046       Standard_Integer nbnodes = GT.NbPoints();
1047       //Check deflection in 2d space for improvement of edge tesselation.
1048       if( isSameParam && nbnodes > 1)
1049       {
1050         Standard_Real aSquareEdgeDef = otherdefedge * otherdefedge;
1051         const TopTools_ListOfShape& lf = theAncestors.FindFromKey(theEdge);
1052         TopTools_ListIteratorOfListOfShape itl(lf);
1053         for (; itl.More(); itl.Next()) {
1054           const TopoDS_Face& aFace = TopoDS::Face (itl.Value());
1055
1056           TopLoc_Location aLoc;
1057           Standard_Real aF, aL;
1058           Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
1059           const Handle(Standard_Type)& aType = aSurf->DynamicType();
1060           if(aType == STANDARD_TYPE(Geom_Plane))
1061             continue;
1062           Handle(Geom2d_Curve) aCurve2d = BRep_Tool::CurveOnSurface(theEdge, aFace, aF, aL);
1063           if(Abs(aF-wFirst)>Precision::PConfusion()||Abs(aL-wLast)>Precision::PConfusion())
1064             continue;
1065
1066           gp_Pnt2d uvf;
1067           Standard_Real parf;
1068           nbnodes = GT.NbPoints();
1069           TColStd_Array1OfReal aParamArray(1, nbnodes);
1070           for (i = 1; i <= nbnodes; i++)
1071           {
1072             GT.Value(cons, theGFace, i, parf, P3d, uvf);
1073             aParamArray.SetValue(i, parf);
1074           }
1075           for (i = 1; i < nbnodes; i++)
1076           {
1077             splitSegment(GT, aSurf, aCurve2d, cons, aSquareEdgeDef, aParamArray(i), aParamArray(i+1), 1);
1078           }
1079         }
1080       }
1081
1082       // Creation of polygons on triangulation:
1083       Standard_Real puv;
1084       nbnodes = GT.NbPoints();
1085
1086       TColStd_Array1OfInteger Nodes(1, nbnodes);
1087       TColStd_Array1OfInteger NodInStruct(1, nbnodes);
1088       TColStd_Array1OfReal Param(1, nbnodes);
1089
1090       // processing of the 1st point:
1091       Nodes(1) = isvf;
1092       NodInStruct(1) = ipf;
1093       Param(1) = wFirst;
1094
1095       // last point:
1096       Nodes(nbnodes) = isvl;
1097       NodInStruct(nbnodes) = ipl;
1098       Param(nbnodes) = wLast;
1099
1100       if(nbnodes > 2)
1101       {
1102         Standard_Integer iv2;
1103         for (i = 2; i < GT.NbPoints(); i++)
1104         {
1105           // Record 3d point
1106           GT.Value(cons, theGFace, i, puv, P3d, uv);
1107           myNbLocat++;
1108           myLocation3d.Bind(myNbLocat, P3d);
1109           NodInStruct(i) = myNbLocat;
1110           // Record 2d point
1111           v2.Initialize(uv.Coord(), myNbLocat, BRepMesh_OnCurve);
1112           iv2=myStructure->AddNode(v2);
1113
1114           Standard_Integer isv = myVemap.FindIndex(iv2);
1115           if (isv == 0) isv = myVemap.Add(iv2);
1116           Nodes(i) = isv;
1117           NodInStruct(i) = myNbLocat;
1118           Param(i) = puv;
1119
1120           if (orEdge == TopAbs_FORWARD)
1121             myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1122           else if (orEdge == TopAbs_REVERSED)
1123             myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1124           else if (orEdge == TopAbs_INTERNAL)
1125             myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1126           ivf = iv2;
1127         }
1128       }
1129
1130       P1 = new Poly_PolygonOnTriangulation(Nodes, Param);
1131       P2 = new Poly_PolygonOnTriangulation(NodInStruct, Param);
1132     }
1133
1134     P2->Deflection(otherdefedge);
1135     BRepMesh_PairOfPolygon pair;
1136     pair.Append(P2);
1137     myEdges.Bind(theEdge, pair);
1138
1139     if (ivf != ivl) {
1140       if (orEdge == TopAbs_FORWARD)
1141         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Frontier));
1142       else if (orEdge == TopAbs_REVERSED)
1143         myStructure->AddLink(BRepMesh_Edge(ivl, ivf, BRepMesh_Frontier));
1144       else if (orEdge == TopAbs_INTERNAL)
1145         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Fixed));
1146     }
1147
1148
1149   }
1150   // If this Edge has been already checked and it is not degenerated,
1151   // the points of the polygon calculated at the first check are retrieved :
1152   else
1153   {
1154     if (degener)
1155     {
1156       // Create 2d polygon for degenerated edge
1157       TColStd_Array1OfInteger Nodes(1, 2);
1158       TColStd_Array1OfReal PPar(1, 2);
1159
1160       Nodes(1) = isvf;
1161       PPar(1) = wFirst;
1162
1163       Nodes(2) = isvl;
1164       PPar(2) = wLast;
1165
1166       P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
1167     }
1168     else
1169     {
1170       // retrieve the polygone:
1171       const BRepMesh_PairOfPolygon& pair = myEdges.Find(theEdge);
1172       const Handle(Poly_PolygonOnTriangulation)& P = pair.First();
1173       const TColStd_Array1OfInteger& NOD = P->Nodes();
1174       Handle(TColStd_HArray1OfReal) Par = P->Parameters();
1175
1176       // creation anew:
1177       const Standard_Integer nbnodes = NOD.Length();
1178       TColStd_Array1OfInteger Nodes(1, nbnodes);
1179       TColStd_Array1OfReal PPar(1, nbnodes);
1180       Standard_Integer iv2;
1181
1182       Nodes(1) = isvf;
1183       PPar(1) = wFirst;
1184
1185       Nodes(nbnodes) = isvl;
1186       PPar(nbnodes) = wLast;
1187
1188       if (nbnodes > 2)
1189       {
1190         Standard_Integer i;
1191         if (BRep_Tool::SameParameter(theEdge))
1192         {
1193           for (i = 2; i < nbnodes; i++)
1194           {
1195             const Standard_Real puv = Par->Value(i);
1196             theC2d->D0(puv, uv);
1197             v2.Initialize(uv.Coord(), NOD(i), BRepMesh_OnCurve);
1198             iv2 = myStructure->AddNode(v2);
1199
1200             Standard_Integer isv = myVemap.FindIndex(iv2);
1201             if (isv == 0) isv = myVemap.Add(iv2);
1202             Nodes(i) = isv;
1203             PPar(i) = puv;
1204
1205             if (orEdge==TopAbs_FORWARD)
1206               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1207             else if (orEdge == TopAbs_REVERSED)
1208               myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1209             else if (orEdge == TopAbs_INTERNAL)
1210               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1211
1212             ivf = iv2;
1213           }
1214         }
1215         else
1216         {
1217           const Standard_Real wFold = Par->Value(Par->Lower());
1218           const Standard_Real wLold = Par->Value(Par->Upper());
1219
1220           Standard_Real wKoef = 1.;
1221           if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
1222           {
1223             wKoef = (wLast - wFirst) / (wLold - wFold);
1224           }
1225
1226           BRepAdaptor_Curve cons(theEdge, theFace);
1227           Extrema_LocateExtPC pcos;
1228           pcos.Initialize(cons, cons.FirstParameter(),
1229               cons.LastParameter(), Precision::PConfusion());
1230
1231           Standard_Real wPrev;
1232           Standard_Real wCur      = wFirst;
1233           Standard_Real wCurFound = wFirst;
1234           for (i = 2; i < nbnodes; i++)
1235           {
1236             P3d = myLocation3d(NOD(i));
1237             // Record 2d point
1238             wPrev = wCur;
1239             wCur  = wFirst + wKoef*(Par->Value(i) - wFold);
1240                   wCurFound += (wCur - wPrev);
1241             pcos.Perform(P3d, wCurFound);
1242             if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
1243             theC2d->D0(wCurFound, uv);
1244             v2.Initialize(uv.Coord(), NOD(i), BRepMesh_OnCurve);
1245             iv2 = myStructure->AddNode(v2);
1246
1247             Standard_Integer isv = myVemap.FindIndex(iv2);
1248             if (isv == 0) isv = myVemap.Add(iv2);
1249                   Nodes(i) = isv;
1250             PPar(i) = wCurFound;
1251
1252             if (orEdge==TopAbs_FORWARD)
1253               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Frontier));
1254             else if (orEdge == TopAbs_REVERSED)
1255               myStructure->AddLink(BRepMesh_Edge(iv2, ivf, BRepMesh_Frontier));
1256             else if (orEdge == TopAbs_INTERNAL)
1257               myStructure->AddLink(BRepMesh_Edge(ivf, iv2, BRepMesh_Fixed));
1258
1259             ivf = iv2;
1260           }
1261         }
1262       }
1263
1264     P1 = new Poly_PolygonOnTriangulation(Nodes, PPar);
1265
1266     if (ivf != ivl) {
1267       if (orEdge == TopAbs_FORWARD)
1268         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Frontier));
1269       else if (orEdge == TopAbs_REVERSED)
1270         myStructure->AddLink(BRepMesh_Edge(ivl, ivf, BRepMesh_Frontier));
1271       else if (orEdge == TopAbs_INTERNAL)
1272         myStructure->AddLink(BRepMesh_Edge(ivf, ivl, BRepMesh_Fixed));
1273       }
1274     }
1275   }
1276
1277   P1->Deflection(theDefEdge);
1278   if (myInternaledges.IsBound(theEdge))
1279   {
1280     BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
1281     if (orEdge == TopAbs_REVERSED)
1282       pair.Append(P1);
1283     else
1284       pair.Prepend(P1);
1285   }
1286   else
1287   {
1288     BRepMesh_PairOfPolygon pair1;
1289     pair1.Append(P1);
1290     myInternaledges.Bind(theEdge, pair1);
1291   }
1292 }
1293
1294
1295 //=======================================================================
1296 //function : Update(edge)
1297 //purpose  :
1298 //=======================================================================
1299 Standard_Boolean BRepMesh_FastDiscret::Update(const TopoDS_Edge&          theEdge,
1300                                               const TopoDS_Face&          theFace,
1301                                               const Handle(Geom2d_Curve)& theC2d,
1302                                               const Standard_Real         theDefEdge,
1303                                               const Standard_Real         theFirst,
1304                                               const Standard_Real         theLast)
1305 {
1306   TopLoc_Location Loc;
1307   Handle(Poly_Triangulation) T, TNull;
1308   Handle(Poly_PolygonOnTriangulation) Poly, NullPoly;
1309
1310   Standard_Integer i = 1;
1311   Standard_Boolean found = Standard_False;
1312   do
1313   {
1314     BRep_Tool::PolygonOnTriangulation(theEdge,Poly,T,Loc,i);
1315     i++;
1316     if (!found && !T.IsNull() && T->HasUVNodes() &&
1317         !Poly.IsNull() && Poly->HasParameters())
1318     {
1319       if (Poly->Deflection() <= 1.1*theDefEdge)
1320       {
1321         // 2d vertex indices
1322         TopAbs_Orientation orEdge = theEdge.Orientation();
1323         Standard_Integer iv1, iv2, ivl;
1324         Standard_Integer isv1, isv, isvl;
1325
1326         // Get range on 2d curve
1327         Standard_Real wFirst, wLast;
1328         BRep_Tool::Range(theEdge, theFace, wFirst, wLast);
1329
1330         // Get end points on 2d curve
1331         gp_Pnt2d uvFirst, uvLast;
1332         BRep_Tool::UVPoints(theEdge, theFace, uvFirst, uvLast);
1333
1334         // Get vertices
1335         TopoDS_Vertex pBegin, pEnd;
1336         TopExp::Vertices(theEdge,pBegin,pEnd);
1337
1338         const Standard_Boolean sameUV =
1339           uvFirst.IsEqual(uvLast, Precision::PConfusion());
1340
1341         //Controle vertice tolerances
1342         BRepAdaptor_Surface  BS(theFace, Standard_False);
1343         Handle(BRepAdaptor_HSurface) gFace = new BRepAdaptor_HSurface(BS);
1344
1345
1346         gp_Pnt pFirst = gFace->Value(uvFirst.X(), uvFirst.Y());
1347         gp_Pnt pLast  = gFace->Value(uvLast.X(), uvLast.Y());
1348
1349         Standard_Real mindist = 10. * Max(pFirst.Distance(BRep_Tool::Pnt(pBegin)),
1350                                           pLast.Distance(BRep_Tool::Pnt(pEnd)));
1351
1352         if (mindist < BRep_Tool::Tolerance(pBegin) ||
1353             mindist < BRep_Tool::Tolerance(pEnd) ) mindist = theDefEdge;
1354
1355         if (sameUV)
1356         {
1357           // 1. is it really sameUV without being degenerated
1358           gp_Pnt2d uvF, uvL;
1359           theC2d->D0(theFirst, uvF);
1360           theC2d->D0(theLast, uvL);
1361           if (!uvFirst.IsEqual(uvF, Precision::PConfusion())) {
1362             uvFirst = uvF;
1363           }
1364           if (!uvLast.IsEqual(uvL, Precision::PConfusion())) {
1365             uvLast = uvL;
1366           }
1367         }
1368
1369         const TColgp_Array1OfPnt&      Nodes   = T->Nodes();
1370         const TColStd_Array1OfInteger& Indices = Poly->Nodes();
1371         Handle(TColStd_HArray1OfReal)  Param   = Poly->Parameters();
1372
1373         const Standard_Integer nbnodes = Indices.Length();
1374         TColStd_Array1OfInteger NewNodes(1, nbnodes);
1375         TColStd_Array1OfInteger NewNodInStruct(1, nbnodes);
1376
1377         gp_Pnt P3d;
1378         gp_XY theUV;
1379
1380         // Process first vertex
1381         Standard_Integer ipf;
1382         if (myVertices.IsBound(pBegin))
1383         {
1384           ipf = myVertices.Find(pBegin);
1385         }
1386         else
1387         {
1388           if (sameUV && myVertices.IsBound(pEnd))
1389           {
1390               ipf = myVertices.Find(pEnd);
1391           }
1392           else
1393           {
1394             P3d = Nodes(Indices(1));
1395             if (!Loc.IsIdentity())
1396               P3d.Transform(Loc.Transformation());
1397             myNbLocat++;
1398             myLocation3d.Bind(myNbLocat,P3d);
1399             ipf = myNbLocat;
1400           }
1401           myVertices.Bind(pBegin,ipf);
1402         }
1403         NewNodInStruct(1) = ipf;
1404
1405         Handle(BRepMesh_FaceAttribute) aFaceAttribute;
1406         GetFaceAttribute ( theFace, aFaceAttribute );
1407         theUV = BRepMesh_FastDiscretFace::FindUV(pBegin, uvFirst, ipf,
1408           gFace, mindist, aFaceAttribute, myLocation2d);
1409
1410         BRepMesh_Vertex vf(theUV,ipf,BRepMesh_Frontier);
1411         iv1 = myStructure->AddNode(vf);
1412         isv1 = myVemap.FindIndex(iv1);
1413         if (isv1 == 0) isv1 = myVemap.Add(iv1);
1414         NewNodes(1) = isv1;
1415
1416         // Process last vertex
1417         Standard_Integer ipl;
1418         if (pEnd.IsSame(pBegin))
1419         {
1420           ipl = ipf;
1421         }
1422         else
1423         {
1424           if (myVertices.IsBound(pEnd))
1425           {
1426             ipl = myVertices.Find(pEnd);
1427           }
1428           else
1429           {
1430             if (sameUV)
1431             {
1432               ipl = ipf;
1433               ivl = iv1;
1434             }
1435             else
1436             {
1437               myNbLocat++;
1438               myLocation3d.Bind(myNbLocat,Nodes(Indices(nbnodes)).Transformed(Loc.Transformation()));
1439               ipl = myNbLocat;
1440             }
1441             myVertices.Bind(pEnd,ipl);
1442           }
1443         }
1444         NewNodInStruct(nbnodes) = ipl;
1445         theUV = BRepMesh_FastDiscretFace::FindUV(pEnd, uvLast, ipl,
1446           gFace, mindist, aFaceAttribute, myLocation2d);
1447
1448         BRepMesh_Vertex vl(theUV,ipl,BRepMesh_Frontier);
1449
1450         ivl = myStructure->AddNode(vl);
1451         isvl = myVemap.FindIndex(ivl);
1452         if (isvl == 0) isvl = myVemap.Add(ivl);
1453
1454         NewNodes(nbnodes) = isvl;
1455
1456         gp_Pnt2d uv;
1457         BRepMesh_Vertex v;
1458
1459         if (BRep_Tool::SameParameter(theEdge))
1460         {
1461           for (i = 2; i < Indices.Length(); i++)
1462           {
1463             // Record 3d point
1464             P3d = Nodes(Indices(i));
1465             if (!Loc.IsIdentity())
1466               P3d.Transform(Loc.Transformation());
1467             myNbLocat++;
1468             myLocation3d.Bind(myNbLocat, P3d);
1469             NewNodInStruct(i) = myNbLocat;
1470             // Record 2d point
1471             uv = theC2d->Value(Param->Value(i));
1472             v.Initialize(uv.Coord(), myNbLocat, BRepMesh_Frontier);
1473             iv2 = myStructure->AddNode(v);
1474             isv = myVemap.FindIndex(iv2);
1475             if (isv == 0) isv = myVemap.Add(iv2);
1476             NewNodes(i) = isv;
1477
1478             //add links
1479             if (orEdge == TopAbs_FORWARD)
1480               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
1481             else if (orEdge == TopAbs_REVERSED)
1482               myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
1483             else if (orEdge == TopAbs_INTERNAL)
1484               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
1485             iv1 = iv2;
1486           }
1487
1488           // last point
1489           if (iv1 != ivl) {
1490             if (orEdge == TopAbs_FORWARD)
1491               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
1492             else if (orEdge == TopAbs_REVERSED)
1493               myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
1494             else if (orEdge == TopAbs_INTERNAL)
1495               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
1496           }
1497
1498
1499         }
1500         else
1501         {
1502           const Standard_Real wFold = Param->Value(Param->Lower());
1503           const Standard_Real wLold = Param->Value(Param->Upper());
1504
1505           Standard_Real wKoef = 1.;
1506           if ((wFold != wFirst || wLold != wLast) && wLold != wFold)
1507           {
1508             wKoef = (wLast - wFirst) / (wLold - wFold);
1509           }
1510
1511           BRepAdaptor_Curve cons(theEdge, theFace);
1512           Extrema_LocateExtPC pcos;
1513           pcos.Initialize(cons, cons.FirstParameter(),
1514                           cons.LastParameter(), Precision::PConfusion());
1515
1516           Standard_Real wPrev;
1517           Standard_Real wCur      = wFirst;
1518           Standard_Real wCurFound = wFirst;
1519           for (i = 2; i < Indices.Length(); i++)
1520           {
1521             // Record 3d point
1522             P3d = Nodes(Indices(i));
1523             if (!Loc.IsIdentity())
1524               P3d.Transform(Loc.Transformation());
1525             myNbLocat++;
1526             myLocation3d.Bind(myNbLocat, P3d);
1527             NewNodInStruct(i) = myNbLocat;
1528             // Record 2d point
1529             wPrev = wCur;
1530             wCur  = wFirst + wKoef*(Param->Value(i) - wFold);
1531             wCurFound += (wCur - wPrev);
1532             pcos.Perform(P3d, wCurFound);
1533             if (pcos.IsDone()) wCurFound = pcos.Point().Parameter();
1534             theC2d->D0(wCurFound, uv);
1535             v.Initialize(uv.Coord(), myNbLocat, BRepMesh_Frontier);
1536             iv2 = myStructure->AddNode(v);
1537             isv = myVemap.FindIndex(iv2);
1538             if (isv == 0) isv = myVemap.Add(iv2);
1539             NewNodes(i) = isv;
1540
1541
1542             //add links
1543             if (orEdge == TopAbs_FORWARD)
1544               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Frontier));
1545             else if (orEdge == TopAbs_REVERSED)
1546               myStructure->AddLink(BRepMesh_Edge(iv2,iv1,BRepMesh_Frontier));
1547             else if (orEdge == TopAbs_INTERNAL)
1548               myStructure->AddLink(BRepMesh_Edge(iv1,iv2,BRepMesh_Fixed));
1549             iv1 = iv2;
1550           }
1551
1552           // last point
1553           if (iv1 != ivl) {
1554             if (orEdge == TopAbs_FORWARD)
1555               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Frontier));
1556             else if (orEdge == TopAbs_REVERSED)
1557               myStructure->AddLink(BRepMesh_Edge(ivl,iv1,BRepMesh_Frontier));
1558             else if (orEdge == TopAbs_INTERNAL)
1559               myStructure->AddLink(BRepMesh_Edge(iv1,ivl,BRepMesh_Fixed));
1560           }
1561         }
1562
1563         Handle(Poly_PolygonOnTriangulation) P1 =
1564           new Poly_PolygonOnTriangulation(NewNodes, Param->Array1());
1565         P1->Deflection(theDefEdge);
1566         if (myInternaledges.IsBound(theEdge))
1567         {
1568           BRepMesh_PairOfPolygon& pair = myInternaledges.ChangeFind(theEdge);
1569           if (theEdge.Orientation() == TopAbs_REVERSED)
1570             pair.Append(P1);
1571           else
1572             pair.Prepend(P1);
1573         }
1574         else
1575         {
1576           BRepMesh_PairOfPolygon pair1;
1577           pair1.Append(P1);
1578           myInternaledges.Bind(theEdge, pair1);
1579         }
1580
1581         Handle(Poly_PolygonOnTriangulation) P2 =
1582           new Poly_PolygonOnTriangulation(NewNodInStruct, Param->Array1());
1583         P2->Deflection(theDefEdge);
1584         BRepMesh_PairOfPolygon pair;
1585         pair.Append(P2);
1586         myEdges.Bind(theEdge, pair);
1587
1588         found = Standard_True;
1589       }
1590       else
1591       {
1592         BRep_Builder B;
1593         B.UpdateEdge(theEdge,NullPoly,T,Loc);
1594         B.UpdateFace(theFace,TNull);
1595       }
1596     }
1597     else if (!T.IsNull() && !T->HasUVNodes())
1598     {
1599       BRep_Builder B;
1600       B.UpdateEdge(theEdge,NullPoly,T,Loc);
1601       B.UpdateFace(theFace,TNull);
1602     }
1603   }
1604   while (!Poly.IsNull());
1605
1606   return found;
1607 }
1608
1609 //=======================================================================
1610 //function : output
1611 //purpose  :
1612 //=======================================================================
1613 Standard_Integer BRepMesh_FastDiscret::NbTriangles() const
1614 {
1615   return myStructure->NbElements();
1616 }
1617
1618 //=======================================================================
1619 //function : Triangle
1620 //purpose  :
1621 //=======================================================================
1622
1623 const BRepMesh_Triangle& BRepMesh_FastDiscret::Triangle
1624   (const Standard_Integer Index) const
1625 {
1626   return myStructure->GetElement(Index);
1627 }
1628
1629 //=======================================================================
1630 //function : NbEdges
1631 //purpose  :
1632 //=======================================================================
1633
1634 Standard_Integer BRepMesh_FastDiscret::NbEdges() const
1635 {
1636   return myStructure->NbLinks();
1637 }
1638
1639 //=======================================================================
1640 //function : Edge
1641 //purpose  :
1642 //=======================================================================
1643
1644 const BRepMesh_Edge& BRepMesh_FastDiscret::Edge(const Standard_Integer Index) const
1645 {
1646   return myStructure->GetLink(Index);
1647 }
1648
1649 //=======================================================================
1650 //function : NbVertices
1651 //purpose  :
1652 //=======================================================================
1653
1654 Standard_Integer BRepMesh_FastDiscret::NbVertices() const
1655 {
1656   return myStructure->NbNodes();
1657 }
1658
1659 //=======================================================================
1660 //function : Vertex
1661 //purpose  :
1662 //=======================================================================
1663
1664 const BRepMesh_Vertex& BRepMesh_FastDiscret::Vertex
1665   (const Standard_Integer Index) const
1666 {
1667   return myStructure->GetNode(Index);
1668 }
1669
1670 //=======================================================================
1671 //function : Pnt
1672 //purpose  :
1673 //=======================================================================
1674
1675 const gp_Pnt& BRepMesh_FastDiscret::Pnt(const Standard_Integer Index) const
1676 {
1677   return myLocation3d(myStructure->GetNode(Index).Location3d());
1678 }
1679
1680 //=======================================================================
1681 //function : VerticesOfDomain
1682 //purpose  :
1683 //=======================================================================
1684
1685 void BRepMesh_FastDiscret::VerticesOfDomain(BRepMesh_MapOfInteger&  Indices) const
1686 {
1687   Indices.Clear();
1688
1689   // recuperate from the map of edges.
1690   const BRepMesh_MapOfInteger& edmap = myStructure->LinkOfDomain();
1691
1692   // iterator on edges.
1693   BRepMesh_MapOfInteger::Iterator iter(edmap);
1694
1695   Standard_Integer ind_edge;
1696   for (iter.Reset(); iter.More(); iter.Next()) {
1697     ind_edge = iter.Key();
1698     const BRepMesh_Edge& Ed = Edge(ind_edge);
1699     Indices.Add(Ed.FirstNode());
1700     Indices.Add(Ed.LastNode());
1701   }
1702 }
1703
1704 BRepMesh_Status BRepMesh_FastDiscret::CurrentFaceStatus() const
1705 {
1706   return myFacestate;
1707 }
1708
1709 //=======================================================================
1710 //function : GetFaceAttribute
1711 //purpose  :
1712 //=======================================================================
1713
1714 Standard_Boolean BRepMesh_FastDiscret::GetFaceAttribute(const TopoDS_Face&              theFace,
1715                                                         Handle(BRepMesh_FaceAttribute)& theFattrib) const
1716 {
1717   if(myMapattrib.IsBound(theFace))
1718   {
1719     theFattrib = myMapattrib(theFace);
1720     return Standard_True;
1721   }
1722   return Standard_False;
1723 }
1724
1725 //=======================================================================
1726 //function : RemoveFaceAttribute
1727 //purpose  :
1728 //=======================================================================
1729
1730 void BRepMesh_FastDiscret::RemoveFaceAttribute(const TopoDS_Face& theFace)
1731 {
1732   if(myMapattrib.IsBound(theFace))
1733     myMapattrib.UnBind(theFace);
1734 }
1735
1736 //=======================================================================
1737 //function : CreateMutexesForSubShapes
1738 //purpose  :
1739 //=======================================================================
1740 void BRepMesh_FastDiscret::CreateMutexesForSubShapes(const TopoDS_Shape& theShape,
1741                                                      const TopAbs_ShapeEnum theType)
1742 {
1743   myMutexProvider.CreateMutexesForSubShapes(theShape, theType);
1744 }
1745
1746 //=======================================================================
1747 //function : RemoveAllMutexes
1748 //purpose  :
1749 //=======================================================================
1750 void BRepMesh_FastDiscret::RemoveAllMutexes()
1751 {
1752   myMutexProvider.RemoveAllMutexes();
1753 }