src/MeshAlgo/MeshAlgo_Delaunay.gxx

   1 // File:        MeshAlgo_Delaunay.gxx
   2 // Created:     Wed May 12 08:58:20 1993
   3 // Author:      Didier PIFFAULT
   4 //              <dpf@zerox>
   5
   6 #include <gp_XY.hxx>
   7 #include <gp_Pnt2d.hxx>
   8 #include <gp_Vec2d.hxx>
   9 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  10 #include <TColStd_ListOfInteger.hxx>
  11 #include <Precision.hxx>
  12 #include <Bnd_Box2d.hxx>
  13 #include <gp.hxx>
  14 #include <TColStd_MapOfInteger.hxx>
  15 #include <MeshDS_MapOfIntegerInteger.hxx>
  16
  17 typedef TColStd_ListIteratorOfListOfInteger  IteratorOnListOfInteger;
  18 typedef TColStd_ListOfInteger                ListOfInteger;
  19
  20 #define EPSEPS Precision::PConfusion()*Precision::PConfusion()
  21
  22 const gp_XY SortingDirection(M_SQRT1_2, M_SQRT1_2);
  23
  24 //#define TRIANGULATION_MESURE
  25
  26 #ifdef TRIANGULATION_MESURE
  27 #include <OSD_Chronometer.hxx>
  28 #include <NCollection_IncAllocator.hxx>
  29 static OSD_Chronometer ChroTrigu;
  30 static OSD_Chronometer ChroSearch;
  31 static OSD_Chronometer ChroDelete;
  32 static OSD_Chronometer ChroDelTri;
  33 static OSD_Chronometer ChroDelCir;
  34 static OSD_Chronometer ChroCreate;
  35 static OSD_Chronometer ChroFront;
  36 Standard_EXPORT Standard_Boolean Triangulation_Chrono;
  37 #endif
  38
  39 //#define TRIANGULATION_DEBUG
  40
  41 #ifdef TRIANGULATION_DEBUG
  42 Standard_IMPORT Standard_Integer Triangulation_Trace;
  43 #endif
  44
  45
  46 //=======================================================================
  47 //function : MeshAlgo_Delaunay
  48 //purpose  :
  49 //=======================================================================
  50 MeshAlgo_Delaunay::MeshAlgo_Delaunay
  51   (MeshAlgo_Array1OfVertex& Vertices, const Standard_Boolean ZPositive)
  52 :  PositiveOrientation(ZPositive), tCircles(Vertices.Length(),new NCollection_IncAllocator())
  53 {
  54   if (Vertices.Length()>2) {
  55     MeshData=new MeshAlgo_DataStructure(new NCollection_IncAllocator(),Vertices.Length());
  56     Init(Vertices);
  57   }
  58 }
  59
  60 //=======================================================================
  61 //function : MeshAlgo_Delaunay
  62 //purpose  :
  63 //=======================================================================
  64 MeshAlgo_Delaunay::MeshAlgo_Delaunay
  65   (const Handle(MeshAlgo_DataStructure)& OldMesh,
  66    MeshAlgo_Array1OfVertex& Vertices,
  67    const Standard_Boolean ZPositive)
  68 :  PositiveOrientation(ZPositive), tCircles(Vertices.Length(),OldMesh->Allocator())
  69 {
  70   MeshData=OldMesh;
  71   if (Vertices.Length()>2)
  72     Init(Vertices);
  73 }
  74
  75
  76 //=======================================================================
  77 //function : Init
  78 //purpose  :
  79 //=======================================================================
  80 void MeshAlgo_Delaunay::Init(MeshAlgo_Array1OfVertex& Vertices)
  81 {
  82   Bnd_Box2d theBox;
  83   TColStd_Array1OfInteger vertexIndices(Vertices.Lower(), Vertices.Upper());
  84   Standard_Integer niver;
  85
  86   for (niver=Vertices.Lower(); niver<=Vertices.Upper(); niver++) {
  87     theBox.Add(gp_Pnt2d(Vertices(niver).Coord()));
  88     vertexIndices(niver)=MeshData->AddNode(Vertices(niver));
  89   }
  90
  91   theBox.Enlarge(Precision::PConfusion());
  92   SuperMesh(theBox);
  93
  94   MeshAlgo_HeapSortIndexedVertex::Sort
  95     (vertexIndices,
  96      MeshAlgo_ComparatorOfIndexedVertex(SortingDirection,
  97                                         Precision::PConfusion(),
  98                                         MeshData));
  99
 100   Compute(vertexIndices);
 101 }
 102
 103
 104 //=======================================================================
 105 //function : MeshAlgo_Delaunay
 106 //purpose  :
 107 //=======================================================================
 108 MeshAlgo_Delaunay::MeshAlgo_Delaunay
 109   (const Handle(MeshAlgo_DataStructure)& OldMesh,
 110    TColStd_Array1OfInteger& VertexIndices,
 111    const Standard_Boolean ZPositive)
 112 :  PositiveOrientation(ZPositive), tCircles(VertexIndices.Length(),OldMesh->Allocator())
 113 {
 114   MeshData=OldMesh;
 115   if (VertexIndices.Length()>2) {
 116     Bnd_Box2d theBox;
 117     Standard_Integer niver;
 118     for (niver=VertexIndices.Lower(); niver<=VertexIndices.Upper(); niver++) {
 119       theBox.Add(gp_Pnt2d(GetVertex(VertexIndices(niver)).Coord()));
 120     }
 121
 122     theBox.Enlarge(Precision::PConfusion());
 123     SuperMesh(theBox);
 124
 125     MeshAlgo_HeapSortIndexedVertex::Sort
 126       (VertexIndices,
 127        MeshAlgo_ComparatorOfIndexedVertex(SortingDirection,
 128                                           Precision::PConfusion(),
 129                                           MeshData));
 130
 131     Compute(VertexIndices);
 132   }
 133 }
 134
 135 //=======================================================================
 136 //function : Compute
 137 //purpose  :
 138 //=======================================================================
 139 void MeshAlgo_Delaunay::Compute (TColStd_Array1OfInteger& VertexIndices)
 140 {
 141   // Insertion of edges of super triangles in the list of free edges:
 142   MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
 143   Standard_Integer e1, e2, e3;
 144   Standard_Boolean o1, o2, o3;
 145   supTrian.Edges(e1, e2, e3, o1, o2, o3);
 146   loopEdges.Bind(e1, Standard_True);
 147   loopEdges.Bind(e2, Standard_True);
 148   loopEdges.Bind(e3, Standard_True);
 149
 150   if (VertexIndices.Length()>0) {
 151
 152     // Creation of 3 triangles with the node and the edges of the super triangle:
 153     Standard_Integer iVert=VertexIndices.Lower();
 154     CreateTriangles(VertexIndices(iVert), loopEdges);
 155
 156     // Insertion of nodes :
 157     Standard_Boolean modif=Standard_True;
 158     Standard_Integer edgeOn, triPerce;
 159
 160     Standard_Integer aVertIdx;
 161     for (iVert++; iVert<=VertexIndices.Upper(); iVert++) {
 162       aVertIdx = VertexIndices(iVert);
 163       const Vertex& refToVert=GetVertex(aVertIdx);
 164       loopEdges.Clear();
 165
 166       // List of indices of circles containing the node :
 167       MeshDS_ListOfInteger& cirL=tCircles.Select(refToVert.Coord());
 168       MeshDS_ListOfInteger::Iterator itT(cirL);
 169
 170       edgeOn=0;
 171       triPerce=0;
 172
 173       for (; itT.More(); itT.Next()) {
 174
 175         // To add a node in the mesh it is necessary to check to conditions
 176         // - the node should be located on the border of the mesh and in an existing triangle
 177         // - all adjacent triangles should belong to a component connected with this triangle
 178         if (Contains(itT.Value(), refToVert, edgeOn)) {
 179           triPerce=itT.Value();
 180           cirL.Remove(itT);
 181           break;
 182         }
 183       }
 184
 185       if (triPerce>0) {
 186         DeleteTriangle(triPerce, loopEdges);
 187
 188         modif=Standard_True;
 189         while (modif && !cirL.IsEmpty()) {
 190           modif=Standard_False;
 191           MeshDS_ListOfInteger::Iterator itT1(cirL);
 192           for (; itT1.More(); itT1.Next()) {
 193             GetTriangle(itT1.Value()).Edges(e1,e2,e3,o1,o2,o3);
 194             if (loopEdges.IsBound(e1) ||
 195                 loopEdges.IsBound(e2) ||
 196                 loopEdges.IsBound(e3)) {
 197               modif=Standard_True;
 198               DeleteTriangle(itT1.Value(), loopEdges);
 199               cirL.Remove(itT1);
 200               break;
 201             }
 202           }
 203         }
 204
 205 #ifdef TRIANGULATION_DEBUG
 206         if (Triangulation_Trace>0) {
 207           cout << " creation of a triangle with the vertex: ";
 208           cout << refToVert.Coord().X() << "  " << refToVert.Coord().Y() << endl;
 209         }
 210 #endif
 211         // Creation of triangles with the current node
 212         // and free edges and removal of these edges from the list of free edges
 213         CreateTriangles(aVertIdx, loopEdges);
 214
 215       }
 216     }
 217
 218     // check that internal edges are not crossed by triangles
 219     MeshDS_MapOfInteger::Iterator itFr(InternalEdges());
 220
 221     // destruction of triancles intersecting internal edges
 222     // and their replacement by makeshift triangles
 223     Standard_Integer nbc;
 224
 225     itFr.Reset();
 226     for (; itFr.More(); itFr.Next()) {
 227       nbc = MeshData->ElemConnectedTo(itFr.Key()).Extent();
 228       if (nbc == 0) {
 229         MeshLeftPolygonOf(itFr.Key(), Standard_True);
 230         MeshLeftPolygonOf(itFr.Key(), Standard_False);
 231       }
 232     }
 233
 234     // adjustment of meshes to boundary edges
 235     FrontierAdjust();
 236
 237   }
 238
 239   // destruction of triangles containing a top of the super triangle
 240   MeshAlgo_SelectorOfDataStructure select(MeshData);
 241   select.NeighboursOfNode(supVert1);
 242   select.NeighboursOfNode(supVert2);
 243   select.NeighboursOfNode(supVert3);
 244   MeshDS_MapOfInteger::Iterator trs(select.Elements());
 245   loopEdges.Clear();
 246   for (;trs.More(); trs.Next()) {
 247     DeleteTriangle(trs.Key(), loopEdges);
 248   }
 249
 250   // all edges that remain free are removed from loopEdges;
 251   // only the boundary edges of the triangulation remain in loopEdges
 252   MeshDS_MapOfIntegerInteger::Iterator itLEd(loopEdges);
 253   for (; itLEd.More(); itLEd.Next()) {
 254     if (MeshData->ElemConnectedTo(itLEd.Key()).IsEmpty())
 255       MeshData->RemoveLink(itLEd.Key());
 256   }
 257
 258   //the tops of the super triangle are destroyed
 259   MeshData->RemoveNode(supVert1);
 260   MeshData->RemoveNode(supVert2);
 261   MeshData->RemoveNode(supVert3);
 262
 263 }
 264
 265 //=======================================================================
 266 //function : ReCompute
 267 //purpose  :
 268 //=======================================================================
 269 void MeshAlgo_Delaunay::ReCompute (TColStd_Array1OfInteger& VertexIndices)
 270 {
 271   MeshData->ClearDomain();
 272
 273   // Initialization of CircleTool :
 274   tCircles.Initialize(VertexIndices.Length());
 275
 276   if (VertexIndices.Length()>2)
 277     Compute(VertexIndices);
 278 }
 279
 280
 281 //=======================================================================
 282 //function : FrontierAdjust
 283 //purpose  :
 284 //=======================================================================
 285 void MeshAlgo_Delaunay::FrontierAdjust()
 286 {
 287   Standard_Integer e1, e2, e3;
 288   Standard_Boolean o1, o2, o3;
 289   MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
 290   MeshDS_ListOfInteger::Iterator itconx;
 291   ListOfInteger tril;
 292
 293   // find external triangles on boundary edges
 294   MeshDS_MapOfInteger::Iterator itFr(Frontier());
 295   for (; itFr.More(); itFr.Next()) {
 296     const MeshDS_PairOfIndex& aPair = MeshData->ElemConnectedTo(itFr.Key());
 297     for(Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; j++) {
 298       const Triangle& trc=GetTriangle(aPair.Index(j));
 299       trc.Edges(e1,e2,e3,o1,o2,o3);
 300       if      ((itFr.Key()==e1 && !o1) ||
 301                (itFr.Key()==e2 && !o2) ||
 302                (itFr.Key()==e3 && !o3))   {
 303 #ifdef TRIANGULATION_DEBUG
 304         if (Triangulation_Trace>0) {
 305           cout << "---> destruction of the triangle " << aPair.Index(j) << endl;
 306         }
 307 #endif
 308         tril.Append(aPair.Index(j));
 309       }
 310     }
 311   }
 312
 313   // destruction  of external triangles on boundary edges
 314   for (; !tril.IsEmpty(); tril.RemoveFirst()) {
 315     DeleteTriangle(tril.First(), loopEdges);
 316   }
 317
 318   // destrucrion of remaining hanging edges :
 319   MeshDS_MapOfIntegerInteger::Iterator itFE(loopEdges);
 320
 321   for (; itFE.More(); itFE.Next()) {
 322     if (MeshData->ElemConnectedTo(itFE.Key()).IsEmpty())
 323       MeshData->RemoveLink(itFE.Key());
 324   }
 325
 326   // destruction of triangles crossing the boundary edges and
 327   // their replacement by makeshift triangles
 328   itFr.Reset();
 329   for (; itFr.More(); itFr.Next()) {
 330     if (MeshData->ElemConnectedTo(itFr.Key()).IsEmpty()) {
 331       MeshLeftPolygonOf(itFr.Key(), Standard_True);
 332     }
 333   }
 334
 335   // After this processing there sometimes remain triangles outside boundaries.
 336   // Destruction of these triangles :
 337   Standard_Integer nbFront;
 338
 339   // For each edge with only one connection
 340   // If one of its tops already passes two boundary edges,
 341   // the connected triangle is outside of the contour
 342   Standard_Boolean again = Standard_True;
 343
 344   while (again) {
 345     tril.Clear();
 346     loopEdges.Clear();
 347
 348     for (itFr.Initialize(FreeEdges()); itFr.More(); itFr.Next()) {
 349       const Edge& edg=GetEdge(itFr.Key());
 350       if (edg.Movability()!=MeshDS_Frontier) {
 351         nbFront=0;
 352         if (MeshData->ElemConnectedTo(itFr.Key()).IsEmpty())
 353           loopEdges.Bind(itFr.Key(), Standard_True);
 354         else {
 355           for (itconx.Init(MeshData->LinkNeighboursOf(edg.FirstNode()));
 356                itconx.More(); itconx.Next()) {
 357             if (GetEdge(itconx.Value()).Movability()==MeshDS_Frontier) {
 358               nbFront++;
 359               if (nbFront>1) break;
 360             }
 361           }
 362           if (nbFront==2) {
 363             const MeshDS_PairOfIndex& aPair = MeshData->ElemConnectedTo(itFr.Key());
 364             for(Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; j++) {
 365               const Standard_Integer elemId = aPair.Index(j);
 366               GetTriangle(elemId).Edges(e1, e2, e3, o1, o2, o3);
 367               if (GetEdge(e1).Movability()==MeshDS_Free &&
 368                   GetEdge(e2).Movability()==MeshDS_Free &&
 369                   GetEdge(e3).Movability()==MeshDS_Free) {
 370 #ifdef TRIANGULATION_DEBUG
 371                 if (Triangulation_Trace>0) {
 372                 cout << " ----> destruction of the triangle" << elemId <<endl;
 373                 }
 374 #endif
 375                 tril.Append(elemId);
 376               }
 377             }
 378           }
 379         }
 380       }
 381     }
 382
 383     // Destruction des triangles :
 384     Standard_Integer kk = 0;
 385     for (; !tril.IsEmpty(); tril.RemoveFirst()) {
 386       DeleteTriangle(tril.First(), loopEdges);
 387       kk++;
 388     }
 389
 390     // destruction of remaining hanging edges
 391     for (itFE.Initialize(loopEdges); itFE.More(); itFE.Next()) {
 392       if (MeshData->ElemConnectedTo(itFE.Key()).IsEmpty())
 393         MeshData->RemoveLink(itFE.Key());
 394     }
 395
 396     if (kk == 0) break;
 397   }
 398
 399   // find external triangles on boundary edges
 400   // because in comlex curved boundaries external triangles can appear
 401   // after "mesh left polygon"
 402    for (itFr.Initialize(Frontier()); itFr.More(); itFr.Next()) {
 403     const MeshDS_PairOfIndex& aPair = MeshData->ElemConnectedTo(itFr.Key());
 404     for(Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; j++) {
 405       const Triangle& trc=GetTriangle(aPair.Index(j));
 406       trc.Edges(e1,e2,e3,o1,o2,o3);
 407       if      ((itFr.Key()==e1 && !o1) ||
 408                (itFr.Key()==e2 && !o2) ||
 409                (itFr.Key()==e3 && !o3))   {
 410 #ifdef TRIANGULATION_DEBUG
 411         if (Triangulation_Trace>0) {
 412           cout << "---> destruction du triangle " << aPair.Index(j) << endl;
 413         }
 414 #endif
 415         tril.Append(aPair.Index(j));
 416       }
 417     }
 418   }
 419
 420   // destruction  of external triangles on boundary edges
 421   for (; !tril.IsEmpty(); tril.RemoveFirst()) {
 422     DeleteTriangle(tril.First(), loopEdges);
 423   }
 424
 425    for (itFE.Initialize(loopEdges); itFE.More(); itFE.Next()) {
 426       if (MeshData->ElemConnectedTo(itFE.Key()).IsEmpty())
 427         MeshData->RemoveLink(itFE.Key());
 428     }
 429
 430
 431   // in some cases there remain unused boundaries check
 432   for (itFr.Initialize(Frontier()); itFr.More(); itFr.Next()) {
 433     if (MeshData->ElemConnectedTo(itFr.Key()).IsEmpty()) {
 434       MeshLeftPolygonOf(itFr.Key(), Standard_True);
 435     }
 436   }
 437 }
 438
 439
 440 //=======================================================================
 441 //function : MeshLeftPolygonOf
 442 //purpose  : Triangulation of the polygon to the left of <indexEdg>.(material side)
 443 //=======================================================================
 444 void MeshAlgo_Delaunay::MeshLeftPolygonOf(const Standard_Integer indexEdg,
 445                                           const Standard_Boolean forwdEdg)
 446 {
 447   const Edge& edg=GetEdge(indexEdg);
 448   TColStd_SequenceOfInteger polyg;
 449   MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
 450   TColStd_MapOfInteger usedEdges;
 451
 452   // Find the polygon
 453   usedEdges.Add(indexEdg);
 454   Standard_Integer debut, prem, pivo;
 455 #ifndef DEB
 456   Standard_Integer ders =0, oth =0;
 457 #else
 458   Standard_Integer ders, oth;
 459 #endif
 460   if (forwdEdg) {
 461     polyg.Append(indexEdg);
 462     prem=edg.FirstNode();
 463     pivo=edg.LastNode();
 464   }
 465   else {
 466     polyg.Append(-indexEdg);
 467     prem=edg.LastNode();
 468     pivo=edg.FirstNode();
 469   }
 470   debut=prem;
 471   const Vertex& debEd=GetVertex(debut);
 472   const Vertex& finEd=GetVertex(pivo);
 473
 474   // Check the presence of the previous edge <indexEdg> :
 475   MeshDS_ListOfInteger::Iterator itLiVer(MeshData->LinkNeighboursOf(prem));
 476   for (; itLiVer.More(); itLiVer.Next()) {
 477     oth=0;
 478     if (itLiVer.Value()!=indexEdg) {
 479       const Edge& nedg=GetEdge(itLiVer.Value());
 480       oth=nedg.LastNode();
 481       if (oth==prem) oth=nedg.FirstNode();
 482       break;
 483     }
 484   }
 485   if (oth==0) {
 486 #ifdef TRIANGULATION_DEBUG
 487     if (Triangulation_Trace>0)
 488       cout << " MeshLeftPolygonOf : No path to the previous Edge!" << endl;
 489 #endif
 490     return;
 491   }
 492
 493
 494   gp_XY vedge(finEd.Coord());
 495   vedge.Subtract(debEd.Coord());
 496   gp_XY ved1(vedge);
 497   gp_XY ved2;
 498   Standard_Integer curEdg=indexEdg, e1, e2, e3;
 499   Standard_Boolean o1, o2, o3;
 500
 501   // Find the nearest to <indexEdg> closed polygon :
 502   Standard_Boolean InMesh, notInters;
 503   Standard_Integer nextedge;
 504   Standard_Real ang, angref;
 505   gp_XY vpivo, vedcur, voth;
 506
 507   while (pivo!=debut) {
 508     nextedge=0;
 509     if (PositiveOrientation) angref=RealFirst();
 510     else                     angref=RealLast();
 511     const Vertex& vertPivo=GetVertex(pivo);
 512     vpivo=vertPivo.Coord();
 513     vpivo.Subtract(debEd.Coord());
 514
 515     itLiVer.Init(MeshData->LinkNeighboursOf(pivo));
 516
 517     // Find the next edge with the greatest angle with <indexEdg>
 518     // and non intersecting <indexEdg> :
 519
 520     for (; itLiVer.More(); itLiVer.Next()) {
 521       if (itLiVer.Value()!=curEdg) {
 522         notInters=Standard_True;
 523         const Edge& nedg=GetEdge(itLiVer.Value());
 524
 525         InMesh=Standard_True;
 526         if (nedg.Movability()==MeshDS_Free) {
 527           if (MeshData->ElemConnectedTo(itLiVer.Value()).IsEmpty())
 528             InMesh=Standard_False;
 529         }
 530
 531         if (InMesh) {
 532           oth=nedg.FirstNode();
 533           if (oth==pivo) oth=nedg.LastNode();
 534
 535           vedcur=GetVertex(oth).Coord();
 536           vedcur.Subtract(vertPivo.Coord());
 537           if (vedcur.Modulus() >= gp::Resolution() &&
 538               ved1.Modulus() >= gp::Resolution()) {
 539             if (prem!=debut && oth!=debut) {
 540               voth=GetVertex(oth).Coord();
 541               voth.Subtract(debEd.Coord());
 542               if ((vedge^vpivo)*(vedge^voth)<0.) {
 543                 voth=vertPivo.Coord();
 544                 voth.Subtract(finEd.Coord());
 545                 if ((vedcur^vpivo)*(vedcur^voth)<0.)
 546                   notInters=Standard_False;
 547               }
 548             }
 549
 550             if (notInters) {
 551               ang=gp_Vec2d(ved1).Angle(gp_Vec2d(vedcur));
 552
 553               if ((PositiveOrientation && ang>angref) ||
 554                   (!PositiveOrientation && ang<angref)) {
 555                 angref=ang;
 556                 ved2=vedcur;
 557                 if (nedg.FirstNode()==pivo) nextedge=itLiVer.Value();
 558                 else                        nextedge=-itLiVer.Value();
 559                 ders=oth;
 560
 561                 //epa: Find correct closing not direct to
 562                 // link but with maximal angle
 563                 // otherwise, polygon greater that expected is found
 564                 //if (ders==debut) break;
 565               }
 566             }
 567           }
 568         }
 569       }
 570     }
 571
 572     if (nextedge!=0) {
 573       if (Abs(nextedge)!=indexEdg && Abs(nextedge)!=curEdg) {
 574         curEdg=Abs(nextedge);
 575
 576         if (!usedEdges.Add(curEdg)) {
 577
 578           //if this edge has already been added to the polygon,
 579           //there is a risk of looping (attention to open contours)
 580 #ifdef TRIANGULATION_DEBUG
 581           if (Triangulation_Trace>0)
 582             cout << " MeshLeftPolygonOf : polygone is not closed!"
 583               << endl;
 584 #endif
 585
 586           // all edges of the boundary contour are removed from the polygon
 587           curEdg=Abs(polyg(polyg.Length()));
 588           while (GetEdge(curEdg).Movability()==MeshDS_Frontier){
 589             polyg.Remove(polyg.Length());
 590             if (polyg.Length()<=0) break;
 591             curEdg=Abs(polyg(polyg.Length()));
 592           }
 593           return;
 594         }
 595
 596         polyg.Append(nextedge);
 597
 598         Standard_Boolean forw=nextedge>0;
 599         const MeshDS_PairOfIndex& aPair = MeshData->ElemConnectedTo(curEdg);
 600         for(Standard_Integer j = 1, jn = aPair.Extent(); j <= jn; j++) {
 601           const Standard_Integer elemId = aPair.Index(j);
 602           GetTriangle(elemId).Edges(e1,e2,e3,o1,o2,o3);
 603           if ((e1==curEdg && o1==forw) ||
 604               (e2==curEdg && o2==forw) ||
 605               (e3==curEdg && o3==forw)) {
 606             DeleteTriangle(elemId, loopEdges);
 607             break;
 608           }
 609         }
 610       }
 611     }
 612     else {
 613 #ifdef TRIANGULATION_DEBUG
 614       if (Triangulation_Trace>0)
 615         cout << " MeshLeftPolygonOf : No next!" << endl;
 616 #endif
 617       return;
 618     }
 619     prem=pivo;
 620     pivo=ders;
 621     ved1=ved2;
 622   }
 623
 624
 625   // Destruction of unused free edges :
 626   MeshDS_MapOfIntegerInteger::Iterator itFE(loopEdges);
 627
 628   for (; itFE.More(); itFE.Next()) {
 629     if (MeshData->ElemConnectedTo(itFE.Key()).IsEmpty())
 630       MeshData->RemoveLink(itFE.Key());
 631   }
 632
 633   MeshPolygon(polyg);
 634 }
 635
 636
 637 //=======================================================================
 638 //function : MeshPolygon
 639 //purpose  : Triangulatiion of a closed polygon described by the list of indexes of
 640 //           its edges in the structure. (negative index == reversed)
 641 //=======================================================================
 642 void MeshAlgo_Delaunay::MeshPolygon  (TColStd_SequenceOfInteger& poly)
 643 {
 644   Standard_Integer vert, vert1, vert2, vert3 =0, tri;
 645
 646   if (poly.Length()==3) {
 647     tri=MeshData->AddElement(Triangle(Abs(poly(1)),Abs(poly(2)),Abs(poly(3)),
 648                                       poly(1)>0,   poly(2)>0,   poly(3)>0,
 649                                       MeshDS_Free));
 650     tCircles.MocAdd(tri);
 651     const Edge& edg1=GetEdge(Abs(poly(1)));
 652     const Edge& edg2=GetEdge(Abs(poly(2)));
 653     if (poly(1)>0) {
 654       vert1=edg1.FirstNode();
 655       vert2=edg1.LastNode();
 656     }
 657     else {
 658       vert1=edg1.LastNode();
 659       vert2=edg1.FirstNode();
 660     }
 661     if (poly(2)>0)
 662       vert3=edg2.LastNode();
 663     else
 664       vert3=edg2.FirstNode();
 665
 666     if (!tCircles.Add(GetVertex(vert1).Coord(),
 667                       GetVertex(vert2).Coord(),
 668                       GetVertex(vert3).Coord(),
 669                       tri)) {
 670       MeshData->RemoveElement(tri);
 671     }
 672   }
 673
 674   else if (poly.Length()>3) {
 675     const Edge& edg=GetEdge(Abs(poly(1)));
 676     Standard_Real distmin=RealLast();
 677     Standard_Integer ip, used =0;
 678
 679     if (poly(1)>0) {
 680       vert1=edg.FirstNode();
 681       vert2=edg.LastNode();
 682     }
 683     else {
 684       vert1=edg.LastNode();
 685       vert2=edg.FirstNode();
 686     }
 687     gp_XY vedg(GetVertex(vert2).Coord()-
 688                GetVertex(vert1).Coord());
 689     Standard_Real modul=vedg.Modulus();
 690     if (modul>0.) {
 691       vedg.SetCoord(vedg.X()/modul, vedg.Y()/modul);
 692
 693       for (ip=3; ip<=poly.Length(); ip++) {
 694         const Edge& nedg=GetEdge(Abs(poly(ip)));
 695         if (poly(ip)>0) vert=nedg.FirstNode();
 696         else            vert=nedg.LastNode();
 697         gp_XY vep(GetVertex(vert).Coord()-GetVertex(vert2).Coord());
 698
 699         Standard_Real dist=vedg^vep;
 700         if (Abs(dist)>Precision::PConfusion()) {
 701           if ((dist>0. && PositiveOrientation) ||
 702               (dist<0. && !PositiveOrientation)) {
 703             if (Abs(dist)<distmin) {
 704               distmin=dist;
 705               vert3=vert;
 706               used=ip;
 707             }
 708           }
 709         }
 710       }
 711     }
 712
 713     Standard_Integer ne2, ne3;
 714     if (distmin<RealLast()) {
 715       ne2=MeshData->AddLink(Edge(vert2, vert3, MeshDS_Free));
 716       ne3=MeshData->AddLink(Edge(vert3, vert1, MeshDS_Free));
 717       tri=MeshData->AddElement(Triangle(Abs(poly(1)), Abs(ne2), Abs(ne3),
 718                                         (poly(1)>0),  (ne2>0),  (ne3>0),
 719                                         MeshDS_Free));
 720
 721       if (!tCircles.Add(GetVertex(vert1).Coord(),
 722                         GetVertex(vert2).Coord(),
 723                         GetVertex(vert3).Coord(),
 724                         tri)) {
 725         MeshData->RemoveElement(tri);
 726       }
 727
 728       if (used<poly.Length()) {
 729         TColStd_SequenceOfInteger suitePoly;
 730         poly.Split(used, suitePoly);
 731         suitePoly.Prepend(-ne3);
 732         MeshPolygon(suitePoly);
 733       }
 734       else
 735         poly.Remove(poly.Length());
 736
 737       if (used>3) {
 738         poly.SetValue(1, -ne2);
 739         MeshPolygon(poly);
 740       }
 741     }
 742     else {
 743 #ifdef TRIANGULATION_DEBUG
 744       if (Triangulation_Trace>0) {
 745         cout << " MeshPolygon : impossible !" << endl;
 746         if (Triangulation_Trace==5) {
 747           cout << " MeshPolygon : ";
 748           for (vert=1; vert<=poly.Length(); vert++)
 749             cout << poly(vert) << " ";
 750           cout<<endl;
 751         }
 752       }
 753 #endif
 754     }
 755   }
 756 }
 757
 758 //=======================================================================
 759 //function : SuperMesh
 760 //purpose  :
 761 //=======================================================================
 762 void MeshAlgo_Delaunay::SuperMesh  (const Bnd_Box2d& theBox)
 763 {
 764   Standard_Real xMin, yMin, xMax, yMax;
 765   theBox.Get(xMin, yMin, xMax, yMax);
 766   Standard_Real deltax=xMax-xMin;
 767   Standard_Real deltay=yMax-yMin;
 768
 769   Standard_Real deltaMin=Min(deltax, deltay);
 770   Standard_Real deltaMax=Max(deltax, deltay);
 771   Standard_Real delta=deltax+deltay;
 772   tCircles.SetMinMaxSize(gp_XY(xMin, yMin), gp_XY(xMax, yMax));
 773
 774   Standard_Integer koeff = 2;
 775   if(MeshData->NbNodes()>100)
 776     koeff = 5;
 777   else if(MeshData->NbNodes()>1000)
 778     koeff = 7;
 779
 780   tCircles.SetCellSize(deltax/koeff, deltay/koeff);
 781
 782   supVert1=MeshData->AddNode(Vertex((xMin+xMax)/2, yMax+deltaMax, MeshDS_Free));
 783   supVert2=MeshData->AddNode(Vertex(xMin-delta, yMin-deltaMin, MeshDS_Free));
 784   supVert3=MeshData->AddNode(Vertex(xMax+delta, yMin-deltaMin, MeshDS_Free));
 785
 786   Standard_Integer niver;
 787   if (!PositiveOrientation) {
 788     niver=supVert2;
 789     supVert2=supVert3;
 790     supVert3=niver;
 791   }
 792
 793   Standard_Integer
 794     ed1=MeshData->AddLink(Edge(supVert1,supVert2,MeshDS_Free));
 795   Standard_Integer
 796     ed2=MeshData->AddLink(Edge(supVert2,supVert3,MeshDS_Free));
 797   Standard_Integer
 798     ed3=MeshData->AddLink(Edge(supVert3,supVert1,MeshDS_Free));
 799   supTrian=Triangle(Abs(ed1), Abs(ed2), Abs(ed3),
 800                     (ed1>0), (ed2>0), (ed3>0), MeshDS_Free);
 801 }
 802
 803
 804 //=======================================================================
 805 //function : CreateTriangles
 806 //purpose  : Creation of triangles from the current node and free edges
 807 //           and deletion of these edges in the list :
 808 //=======================================================================
 809 void MeshAlgo_Delaunay::CreateTriangles (const Standard_Integer theVertexIndex,
 810                                          MeshDS_MapOfIntegerInteger& theEdges)
 811 {
 812   MeshDS_MapOfIntegerInteger::Iterator itFE(theEdges);
 813   Standard_Real z12, z23, modul;
 814   Standard_Integer ne1, ne3, nt;
 815   gp_XY vl1, vl2, vl3;
 816   ListOfInteger EdgLoop, EdgOK, EdgExter;
 817
 818   const gp_XY& VertexCoord = MeshData->GetNode(theVertexIndex).Coord();
 819   for (; itFE.More(); itFE.Next()) {
 820     const Edge& edg=GetEdge(itFE.Key());
 821     Standard_Integer deb=edg.FirstNode();
 822     Standard_Integer fin=edg.LastNode();
 823     Standard_Boolean sens=(Standard_Boolean)theEdges(itFE.Key());
 824     if (!sens) {
 825       nt=deb;
 826       deb=fin;
 827       fin=nt;
 828     }
 829
 830     const Vertex& debVert=GetVertex(deb);
 831     const Vertex& finVert=GetVertex(fin);
 832
 833     vl1=debVert.Coord();
 834     vl1.Subtract(VertexCoord);
 835     vl2=finVert.Coord();
 836     vl2.Subtract(debVert.Coord());
 837 //    vl3=theVertex.Coord();
 838     vl3=VertexCoord;
 839     vl3.Subtract(finVert.Coord());
 840
 841     z12=z23=0.;
 842     modul=vl2.Modulus();
 843     if (modul>Precision::PConfusion()) {
 844       vl2.SetCoord(vl2.X()/modul, vl2.Y()/modul);
 845       z12=vl1^vl2;
 846       z23=vl2^vl3;
 847     }
 848
 849     if (Abs(z12)>=Precision::PConfusion() &&
 850         Abs(z23)>=Precision::PConfusion()) {
 851       Standard_Boolean sensOK;
 852       if (PositiveOrientation) sensOK=(z12>0. && z23>0.);
 853       else                     sensOK=(z12<0. && z23<0.);
 854       if (sensOK) {
 855         ne1=MeshData->AddLink(Edge(theVertexIndex, deb, MeshDS_Free));
 856         ne3=MeshData->AddLink(Edge(fin, theVertexIndex, MeshDS_Free));
 857         nt=MeshData->AddElement(Triangle(Abs(ne1), itFE.Key(), Abs(ne3),
 858                                          (ne1>0), sens, (ne3>0),
 859                                          MeshDS_Free));
 860
 861         if (!tCircles.Add(VertexCoord,
 862                           debVert.Coord(),
 863                           finVert.Coord(), nt))
 864           MeshData->RemoveElement(nt);
 865       }
 866       else {
 867
 868         if (sens) EdgLoop.Append(itFE.Key());
 869         else      EdgLoop.Append(-itFE.Key());
 870         if (vl1.SquareModulus()>vl3.SquareModulus()) {
 871           ne1=MeshData->AddLink(Edge(theVertexIndex, deb, MeshDS_Free));
 872           EdgExter.Append(Abs(ne1));
 873         }
 874         else{
 875           ne3=MeshData->AddLink(Edge(fin, theVertexIndex, MeshDS_Free));
 876           EdgExter.Append(Abs(ne3));
 877         }
 878       }
 879     }
 880 #ifdef TRIANGULATION_DEBUG
 881     else {
 882       if (Triangulation_Trace>0)
 883         cout << " CreateTriangles : too small vector product !" << endl;
 884     }
 885 #endif
 886   }
 887   theEdges.Clear();
 888   while (!EdgExter.IsEmpty()) {
 889     const MeshDS_PairOfIndex& conx = MeshData->ElemConnectedTo(Abs(EdgExter.First()));
 890     if (!conx.IsEmpty())
 891       DeleteTriangle(conx.FirstIndex(), theEdges);
 892     EdgExter.RemoveFirst();
 893   }
 894
 895   for (itFE.Initialize(theEdges); itFE.More(); itFE.Next()) {
 896     if (MeshData->ElemConnectedTo(itFE.Key()).IsEmpty())
 897       MeshData->RemoveLink(itFE.Key());
 898   }
 899
 900   while (!EdgLoop.IsEmpty()) {
 901     if (GetEdge(Abs(EdgLoop.First())).Movability()!=MeshDS_Deleted) {
 902       MeshLeftPolygonOf(Abs(EdgLoop.First()), (EdgLoop.First()>0));
 903     }
 904     EdgLoop.RemoveFirst();
 905   }
 906 }
 907
 908 //=======================================================================
 909 //function : DeleteTriangle
 910 //purpose : The concerned triangles are deleted and the freed edges are added in
 911 //           <loopEdges>.  If an edge is added twice, it does not exist and
 912 //          it is necessary to destroy it.  This corresponds to the destruction of two
 913 //          connected triangles.
 914 //=======================================================================
 915
 916 void  MeshAlgo_Delaunay::DeleteTriangle (const Standard_Integer tIndex,
 917                                          MeshDS_MapOfIntegerInteger& fEdges)
 918 {
 919   tCircles.Delete(tIndex);
 920
 921   Standard_Integer fe1, fe2, fe3;
 922   Standard_Boolean or1, or2, or3;
 923   GetTriangle(tIndex).Edges(fe1, fe2, fe3, or1, or2, or3);
 924   MeshData->RemoveElement(tIndex);
 925
 926   if (!fEdges.Bind(fe1, or1)) {
 927     fEdges.UnBind(fe1);
 928     MeshData->RemoveLink(fe1);
 929   }
 930   if (!fEdges.Bind(fe2, or2))  {
 931     fEdges.UnBind(fe2);
 932     MeshData->RemoveLink(fe2);
 933   }
 934   if (!fEdges.Bind(fe3, or3))  {
 935     fEdges.UnBind(fe3);
 936     MeshData->RemoveLink(fe3);
 937   }
 938
 939 }
 940
 941
 942 //=======================================================================
 943 //function : AddVertex
 944 //purpose  :
 945 //=======================================================================
 946 void  MeshAlgo_Delaunay::AddVertex(const Vertex& theVert)
 947 {
 948   Standard_Integer nv = MeshData->AddNode(theVert);
 949
 950   // Iterator in the list of indexes of circles containing the node :
 951   MeshDS_ListOfInteger& cirL=tCircles.Select(theVert.Coord());
 952
 953   Standard_Integer edgon=0;
 954   Standard_Integer triPer=0;
 955   Standard_Integer e1, e2, e3;
 956   Standard_Boolean o1, o2, o3;
 957
 958   MeshDS_ListOfInteger::Iterator itT(cirL);
 959   for (; itT.More(); itT.Next()) {
 960
 961     // To add a node in the mesh it is necessary to check conditions:
 962     // - the node should be within the boundaries of the mesh and so in an existing triangle
 963     // - all adjacent triangles should belong to a component connected with this triangle
 964     if (Contains(itT.Value(), theVert, edgon)) {
 965       if (edgon==0) {
 966         triPer=itT.Value();
 967         cirL.Remove(itT);
 968         break;
 969       }
 970       else if (GetEdge(edgon).Movability()==MeshDS_Free) {
 971         triPer=itT.Value();
 972         cirL.Remove(itT);
 973         break;
 974       }
 975     }
 976   }
 977
 978   if (triPer>0) {
 979
 980     MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
 981     DeleteTriangle(triPer, loopEdges);
 982
 983     Standard_Boolean modif=Standard_True;
 984     while (modif && !cirL.IsEmpty()) {
 985       modif=Standard_False;
 986       MeshDS_ListOfInteger::Iterator itT1(cirL);
 987       for (; itT1.More(); itT1.Next()) {
 988         GetTriangle(itT.Value()).Edges(e1,e2,e3,o1,o2,o3);
 989         if (loopEdges.IsBound(e1) ||
 990             loopEdges.IsBound(e2) ||
 991             loopEdges.IsBound(e3)) {
 992           modif=Standard_True;
 993           DeleteTriangle(itT1.Value(), loopEdges);
 994           cirL.Remove(itT1);
 995           break;
 996         }
 997       }
 998     }
 999
1000     // Creation of triangles with the current node and free edges
1001     // and removal of these edges from the list of free edges
1002     CreateTriangles(nv, loopEdges);
1003
1004     // Check that internal edges are not crossed by the triangles
1005     MeshDS_MapOfInteger::Iterator itFr(InternalEdges());
1006
1007     // Destruction of triangles crossing internal edges and
1008     // their replacement by makeshift triangles
1009     Standard_Integer nbc;
1010     itFr.Reset();
1011     for (; itFr.More(); itFr.Next()) {
1012       nbc = MeshData->ElemConnectedTo(itFr.Key()).Extent();
1013       if (nbc == 0) {
1014         MeshLeftPolygonOf(itFr.Key(), Standard_True);
1015         MeshLeftPolygonOf(itFr.Key(), Standard_False);
1016       }
1017     }
1018
1019     FrontierAdjust();
1020
1021   }
1022
1023 }
1024
1025 //=======================================================================
1026 //function : RemoveVertex
1027 //purpose  :
1028 //=======================================================================
1029 void  MeshAlgo_Delaunay::RemoveVertex(const Vertex& theVert)
1030 {
1031   MeshAlgo_SelectorOfDataStructure select(MeshData);
1032   select.NeighboursOf(theVert);
1033
1034   MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
1035
1036   // Loop on triangles to be destroyed :
1037   MeshDS_MapOfInteger::Iterator trs(select.Elements());
1038   for (;trs.More(); trs.Next()) {
1039     DeleteTriangle(trs.Key(), loopEdges);
1040   }
1041
1042   TColStd_SequenceOfInteger polyg;
1043   Standard_Integer iv;
1044   Standard_Integer nbLi=loopEdges.Extent();
1045   MeshDS_MapOfIntegerInteger::Iterator itFE(loopEdges);
1046
1047   if (itFE.More()) {
1048     const Edge& edg=GetEdge(itFE.Key());
1049     Standard_Integer deb=edg.FirstNode();
1050     Standard_Integer fin;
1051     Standard_Integer pivo=edg.LastNode();
1052     Standard_Integer iseg=itFE.Key();
1053     Standard_Boolean sens=(Standard_Boolean)loopEdges(iseg);
1054     if (!sens) {
1055       iv=deb;
1056       deb=pivo;
1057       pivo=iv;
1058       polyg.Append(-iseg);
1059     }
1060     else {
1061       polyg.Append(iseg);
1062     }
1063     loopEdges.UnBind(iseg);
1064     fin=deb;
1065     MeshDS_ListOfInteger::Iterator itLiV;
1066     Standard_Integer vcur;
1067     while (pivo!=fin) {
1068       itLiV.Init(MeshData->LinkNeighboursOf(pivo));
1069       for (; itLiV.More(); itLiV.Next()) {
1070         if (itLiV.Value()!=iseg && loopEdges.IsBound(itLiV.Value())) {
1071           iseg=itLiV.Value();
1072           const Edge& edg1=GetEdge(iseg);
1073           vcur=edg1.LastNode();
1074           if (vcur!=pivo) {
1075             vcur=edg1.FirstNode();
1076             polyg.Append(-iseg);
1077           }
1078           else
1079             polyg.Append(iseg);
1080           pivo=vcur;
1081           loopEdges.UnBind(iseg);
1082           break;
1083         }
1084       }
1085       if (nbLi<=0) break;
1086       nbLi--;
1087     }
1088     MeshPolygon(polyg);
1089   }
1090 }
1091
1092
1093 //=======================================================================
1094 //function : AddVertices
1095 //purpose  :
1096 //=======================================================================
1097 void  MeshAlgo_Delaunay::AddVertices(MeshAlgo_Array1OfVertex& vertices)
1098 {
1099   MeshAlgo_HeapSortVertex::Sort
1100     (vertices,
1101      MeshAlgo_ComparatorOfVertex(SortingDirection, Precision::PConfusion()));
1102
1103   MeshDS_MapOfIntegerInteger loopEdges(10,MeshData->Allocator());
1104   Standard_Boolean modif=Standard_True;
1105   Standard_Integer edgon, triPer;
1106   Standard_Integer e1, e2, e3;
1107   Standard_Boolean o1, o2, o3;
1108
1109   Standard_Integer niver;
1110   Standard_Integer aIdxVert;
1111   for (niver=vertices.Lower(); niver<=vertices.Upper(); niver++) {
1112     aIdxVert = MeshData->AddNode(vertices(niver));
1113
1114     // Iterator in the list of indexes of circles containing the node
1115     MeshDS_ListOfInteger& cirL=tCircles.Select(vertices(niver).Coord());
1116
1117     edgon=0;
1118     triPer=0;
1119
1120     MeshDS_ListOfInteger::Iterator itT(cirL);
1121     for (; itT.More(); itT.Next()) {
1122
1123       // To add a node in the mesh it is necessary to check conditions:
1124       // the node should be within the boundaries of the mesh and so in an existing triangle
1125       // all adjacent triangles should belong to a component connected with this triangle
1126       if (Contains(itT.Value(), vertices(niver), edgon)) {
1127         if (edgon==0) {
1128           triPer=itT.Value();
1129           cirL.Remove(itT);
1130           break;
1131         }
1132         else if (GetEdge(edgon).Movability()==MeshDS_Free) {
1133           triPer=itT.Value();
1134           cirL.Remove(itT);
1135           break;
1136         }
1137       }
1138     }
1139
1140     if (triPer>0) {
1141       DeleteTriangle(triPer, loopEdges);
1142
1143       modif=Standard_True;
1144       while (modif && !cirL.IsEmpty()) {
1145         modif=Standard_False;
1146         MeshDS_ListOfInteger::Iterator itT1(cirL);
1147         for (; itT1.More(); itT1.Next()) {
1148           GetTriangle(itT1.Value()).Edges(e1,e2,e3,o1,o2,o3);
1149           if (loopEdges.IsBound(e1) ||
1150               loopEdges.IsBound(e2) ||
1151               loopEdges.IsBound(e3)) {
1152             modif=Standard_True;
1153             DeleteTriangle(itT1.Value(), loopEdges);
1154             cirL.Remove(itT1);
1155             break;
1156           }
1157         }
1158       }
1159
1160       // Creation of triangles with the current node and free edges
1161       // and removal of these edges from the list of free edges
1162       CreateTriangles(aIdxVert, loopEdges);
1163     }
1164   }
1165
1166     // Check that internal edges are not crossed by triangles
1167     MeshDS_MapOfInteger::Iterator itFr(InternalEdges());
1168
1169     // Destruction of triangles crossing internal edges
1170     //and their replacement by makeshift triangles
1171     Standard_Integer nbc;
1172     itFr.Reset();
1173     for (; itFr.More(); itFr.Next()) {
1174       nbc = MeshData->ElemConnectedTo(itFr.Key()).Extent();
1175       if (nbc == 0) {
1176         MeshLeftPolygonOf(itFr.Key(), Standard_True);
1177         MeshLeftPolygonOf(itFr.Key(), Standard_False);
1178       }
1179     }
1180
1181   // Adjustment of meshes to boundary edges
1182   FrontierAdjust();
1183 }
1184
1185
1186 //=======================================================================
1187 //function : RevertDiagonal
1188 //purpose  :
1189 //=======================================================================
1190 Standard_Boolean MeshAlgo_Delaunay::RevertDiagonal(const Standard_Integer ind)
1191 {
1192   const MeshDS_PairOfIndex& elConx = MeshData->ElemConnectedTo(ind);
1193   const Edge& lEdge = GetEdge(ind);
1194   if (elConx.Extent()==2 && lEdge.Movability()==MeshDS_Free) {
1195     Standard_Integer t1(elConx.FirstIndex());
1196     Standard_Integer t2(elConx.LastIndex());
1197
1198     Standard_Integer e1t1, e2t1, e3t1, e1t2, e2t2, e3t2 ;
1199     Standard_Boolean o1t1, o2t1, o3t1, o1t2, o2t2, o3t2;
1200 #ifndef DEB
1201     Standard_Integer ed13=0, ed23=0, ed14=0, ed24=0, v1, v2, v3=0, v4=0, vc1;
1202     Standard_Boolean oindt1=Standard_False, or13=Standard_False,
1203     or23=Standard_False, or14=Standard_False, or24=Standard_False, orien;
1204 #else
1205     Standard_Integer ed13, ed23, ed14, ed24, v1, v2, v3, v4, vc1;
1206     Standard_Boolean oindt1, or13, or23, or14, or24, orien;
1207 #endif
1208     GetTriangle(t1).Edges(e1t1, e2t1, e3t1, o1t1, o2t1, o3t1);
1209     GetTriangle(t2).Edges(e1t2, e2t2, e3t2, o1t2, o2t2, o3t2);
1210
1211     v1=lEdge.FirstNode(); v2=lEdge.LastNode();
1212     if      (e1t1==ind) {
1213       if (o2t1) v3 =GetEdge(e2t1).LastNode();
1214       else      v3 =GetEdge(e2t1).FirstNode();
1215       ed13=e3t1; ed23=e2t1;
1216       or13=o3t1; or23=o2t1;
1217       oindt1=o1t1;
1218     }
1219     else if (e2t1==ind) {
1220       if (o3t1) v3 =GetEdge(e3t1).LastNode();
1221       else      v3 =GetEdge(e3t1).FirstNode();
1222       ed13=e1t1; ed23=e3t1;
1223       or13=o1t1; or23=o3t1;
1224       oindt1=o2t1;
1225     }
1226     else if (e3t1==ind) {
1227       if (o1t1) v3 =GetEdge(e1t1).LastNode();
1228       else      v3 =GetEdge(e1t1).FirstNode();
1229       ed13=e2t1; ed23=e1t1;
1230       or13=o2t1; or23=o1t1;
1231       oindt1=o3t1;
1232     }
1233     if      (e1t2==ind) {
1234       if (o2t2) v4 =GetEdge(e2t2).LastNode();
1235       else      v4 =GetEdge(e2t2).FirstNode();
1236       ed14=e2t2; ed24=e3t2;
1237       or14=o2t2; or24=o3t2;
1238     }
1239     else if (e2t2==ind) {
1240       if (o3t2) v4 =GetEdge(e3t2).LastNode();
1241       else      v4 =GetEdge(e3t2).FirstNode();
1242       ed14=e3t2; ed24=e1t2;
1243       or14=o3t2; or24=o1t2;
1244     }
1245     else if (e3t2==ind) {
1246       if (o1t2) v4 =GetEdge(e1t2).LastNode();
1247       else      v4 =GetEdge(e1t2).FirstNode();
1248       ed14=e1t2; ed24=e2t2;
1249       or14=o1t2; or24=o2t2;
1250     }
1251     if (!oindt1) {
1252       vc1=v3; v3=v4; v4=vc1;
1253       vc1=ed13; ed13=ed24; ed24=vc1;
1254       orien =or13; or13=or24; or24=orien ;
1255       vc1=ed14; ed14=ed23; ed23=vc1;
1256       orien =or14; or14=or23; or23=orien ;
1257     }
1258     const Vertex& vert1 = GetVertex(v1);
1259     const Vertex& vert2 = GetVertex(v2);
1260     const Vertex& vert3 = GetVertex(v3);
1261     const Vertex& vert4 = GetVertex(v4);
1262
1263     gp_XY ved13(vert1.Coord()); ved13.Subtract(vert3.Coord());
1264     gp_XY ved14(vert4.Coord()); ved14.Subtract(vert1.Coord());
1265     gp_XY ved23(vert3.Coord()); ved23.Subtract(vert2.Coord());
1266     gp_XY ved24(vert2.Coord()); ved24.Subtract(vert4.Coord());
1267
1268     Standard_Real z13, z24, modul;
1269     z13=z24=0.;
1270     modul=ved13.Modulus();
1271     if (modul>Precision::PConfusion()) {
1272       ved13.SetCoord(ved13.X()/modul, ved13.Y()/modul);
1273       z13=ved13^ved14;
1274     }
1275     modul=ved24.Modulus();
1276     if (modul>Precision::PConfusion()) {
1277       ved24.SetCoord(ved24.X()/modul, ved24.Y()/modul);
1278       z24=ved24^ved23;
1279     }
1280
1281     if (Abs(z13)>=Precision::PConfusion()&&Abs(z24)>=Precision::PConfusion()) {
1282       if ((z13>0. && z24>0.) || (z13<0. && z24<0.)) {
1283         tCircles.Delete(t1);
1284         tCircles.Delete(t2);
1285         if (!tCircles.Add(vert4.Coord(), vert2.Coord(), vert3.Coord(), t1) &&
1286             !tCircles.Add(vert3.Coord(), vert1.Coord(), vert4.Coord(), t2)) {
1287           Standard_Integer newd=ind;
1288           Edge newEdg=Edge(v3, v4, MeshDS_Free);
1289           if (!MeshData->SubstituteLink(newd, newEdg)) {
1290             newd=MeshData->IndexOf(newEdg);
1291             MeshData->RemoveLink(ind);
1292           }
1293           MeshData->SubstituteElement(t1, Triangle(ed24, ed23, newd,
1294                                                    or24, or23, Standard_True,
1295                                                    MeshDS_Free));
1296           MeshData->SubstituteElement(t2, Triangle(ed13, ed14, newd,
1297                                                    or13, or14, Standard_False,
1298                                                    MeshDS_Free));
1299           return Standard_True;
1300         }
1301         else {
1302           if (oindt1) {
1303             tCircles.Add(vert1.Coord(), vert2.Coord(), vert3.Coord(), t1);
1304             tCircles.Add(vert2.Coord(), vert1.Coord(), vert4.Coord(), t2);
1305           }
1306           else {
1307             tCircles.Add(vert1.Coord(), vert2.Coord(), vert3.Coord(), t2);
1308             tCircles.Add(vert2.Coord(), vert1.Coord(), vert4.Coord(), t1);
1309           }
1310         }
1311       }
1312     }
1313   }
1314   return Standard_False;
1315 }
1316
1317 //=======================================================================
1318 //function : UseEdge
1319 //purpose  :
1320 //=======================================================================
1321 Standard_Boolean MeshAlgo_Delaunay::UseEdge(const Standard_Integer ind)
1322 {
1323   const MeshDS_PairOfIndex& elConx=MeshData->ElemConnectedTo(ind);
1324   if (elConx.Extent()==0) {
1325     const Edge& lEdge = GetEdge(ind);
1326     Standard_Integer vdeb, pivo, othV, leftEdge, rightEdge;
1327     vdeb=lEdge.FirstNode();
1328     pivo=lEdge.LastNode();
1329     const MeshDS_ListOfInteger& neigVDeb = MeshData->LinkNeighboursOf(vdeb);
1330     const MeshDS_ListOfInteger& neigPivo = MeshData->LinkNeighboursOf(pivo);
1331     if (neigVDeb.Extent()>0 && neigPivo.Extent()>0) {
1332       const Vertex& vertDeb=GetVertex(vdeb);
1333       const Vertex& vertPivo=GetVertex(pivo);
1334
1335       gp_XY vedcur;
1336       gp_XY vedge(vertPivo.Coord());
1337       vedge.Subtract(vertDeb.Coord());
1338
1339       MeshDS_ListOfInteger::Iterator itNeig(neigPivo);
1340 #ifndef DEB
1341       Standard_Real ang =0.;
1342 #else
1343       Standard_Real ang;
1344 #endif
1345       Standard_Real angMin=RealLast();
1346       Standard_Real angMax=RealFirst();
1347       Standard_Boolean InMesh;
1348       leftEdge=rightEdge=0;
1349
1350       for (; itNeig.More(); itNeig.Next()) {
1351         if (itNeig.Value()!=ind) {
1352           const Edge& nedg=GetEdge(itNeig.Value());
1353
1354           InMesh=Standard_True;
1355           if (nedg.Movability()==MeshDS_Free) {
1356             if (MeshData->ElemConnectedTo(itNeig.Value()).IsEmpty())
1357               InMesh=Standard_False;
1358           }
1359
1360           if (InMesh) {
1361             othV=nedg.FirstNode();
1362             if (othV==pivo) othV=nedg.LastNode();
1363
1364             vedcur=GetVertex(othV).Coord();
1365             vedcur.Subtract(vertPivo.Coord());
1366
1367             ang=gp_Vec2d(vedge).Angle(gp_Vec2d(vedcur));
1368           }
1369           if (ang>angMax) {
1370             angMax=ang;
1371             leftEdge=itNeig.Value();
1372           }
1373           if (ang<angMin) {
1374             angMin=ang;
1375             rightEdge=itNeig.Value();
1376           }
1377         }
1378       }
1379       if (leftEdge>0) {
1380         if (leftEdge==rightEdge) {
1381         }
1382         else {
1383         }
1384       }
1385     }
1386   }
1387   return Standard_False;
1388 }
1389
1390 //=======================================================================
1391 //function : SmoothMesh
1392 //purpose  :
1393 //=======================================================================
1394 void MeshAlgo_Delaunay::SmoothMesh(const Standard_Real Epsilon)
1395 {
1396   Standard_Integer baryVert, polyVert, nbPolyVert;
1397   Standard_Real uSom, vSom, newU, newV;
1398   Standard_Integer nbVert=MeshData->NbNodes();
1399   MeshDS_ListOfInteger::Iterator itNeig;
1400
1401   uSom=vSom=0;
1402   for (baryVert=1; baryVert<=nbVert; baryVert++) {
1403     const Vertex& curVert=GetVertex(baryVert);
1404     if (curVert.Movability()==MeshDS_Free) {
1405       const MeshDS_ListOfInteger& neighEdg=MeshData->LinkNeighboursOf(baryVert);
1406       if (neighEdg.Extent()>2) {
1407         nbPolyVert=0;
1408         for (itNeig.Init(neighEdg); itNeig.More(); itNeig.Next()) {
1409           const Edge& nedg=GetEdge(itNeig.Value());
1410           polyVert=nedg.FirstNode();
1411           if (polyVert==baryVert) polyVert=nedg.LastNode();
1412           nbPolyVert++;
1413           const gp_XY& pVal = GetVertex(polyVert).Coord();
1414           uSom+=pVal.X();
1415           vSom+=pVal.Y();
1416         }
1417         if (nbPolyVert>2) {
1418           newU=uSom/(Standard_Real)nbPolyVert;
1419           newV=vSom/(Standard_Real)nbPolyVert;
1420           if (!curVert.Coord().IsEqual(gp_XY(newU, newV), Epsilon)) {
1421             Vertex newVert(newU, newV, curVert.Movability());
1422             MeshData->MoveNode(baryVert, newVert);
1423           }
1424         }
1425       }
1426     }
1427   }
1428 }
1429
1430 //=======================================================================
1431 //function : Result
1432 //purpose  :
1433 //=======================================================================
1434 const Handle(MeshAlgo_DataStructure)& MeshAlgo_Delaunay::Result()const
1435 {
1436   return MeshData;
1437 }
1438
1439 //=======================================================================
1440 //function : Frontier
1441 //purpose  :
1442 //=======================================================================
1443 const MeshDS_MapOfInteger& MeshAlgo_Delaunay::Frontier ()
1444 {
1445   MeshDS_MapOfInteger::Iterator triDom(MeshData->LinkOfDomain());
1446
1447   mapEdges.Clear();
1448   for (; triDom.More(); triDom.Next()) {
1449     if (GetEdge(triDom.Key()).Movability()==MeshDS_Frontier) {
1450       mapEdges.Add(triDom.Key());
1451     }
1452   }
1453   return mapEdges;
1454 }
1455
1456 //=======================================================================
1457 //function : InternalEdges
1458 //purpose  :
1459 //=======================================================================
1460 const MeshDS_MapOfInteger& MeshAlgo_Delaunay::InternalEdges ()
1461 {
1462   MeshDS_MapOfInteger::Iterator triDom(MeshData->LinkOfDomain());
1463
1464   mapEdges.Clear();
1465   for (; triDom.More(); triDom.Next()) {
1466     if (GetEdge(triDom.Key()).Movability()==MeshDS_Fixed) {
1467       mapEdges.Add(triDom.Key());
1468     }
1469   }
1470   return mapEdges;
1471 }
1472
1473 //=======================================================================
1474 //function : FreeEdges
1475 //purpose  :
1476 //=======================================================================
1477 const MeshDS_MapOfInteger& MeshAlgo_Delaunay::FreeEdges ()
1478 {
1479   MeshDS_MapOfInteger::Iterator triDom(MeshData->LinkOfDomain());
1480
1481   mapEdges.Clear();
1482   for (; triDom.More(); triDom.Next()) {
1483     if (MeshData->ElemConnectedTo(triDom.Key()).Extent()<=1) {
1484       mapEdges.Add(triDom.Key());
1485     }
1486   }
1487   return mapEdges;
1488 }
1489
1490
1491 //=======================================================================
1492 //function : Contains
1493 //purpose  :
1494 //=======================================================================
1495 Standard_Boolean MeshAlgo_Delaunay::Contains(const Standard_Integer tri,
1496                                                   const Vertex& vert,
1497                                                   Standard_Integer& edgOn)const
1498 {
1499   edgOn=0;
1500   Standard_Integer e1, e2, e3, p1, p2, p3;
1501   Standard_Boolean o1, o2, o3;
1502   GetTriangle(tri).Edges(e1, e2, e3, o1, o2, o3);
1503   const Edge& edg1=GetEdge(e1);
1504   const Edge& edg2=GetEdge(e2);
1505   const Edge& edg3=GetEdge(e3);
1506   if (o1) {
1507     p1=edg1.FirstNode();
1508     p2=edg1.LastNode();
1509   }
1510   else {
1511     p2=edg1.FirstNode();
1512     p1=edg1.LastNode();
1513   }
1514   if (o3) p3=edg3.FirstNode();
1515   else    p3=edg3.LastNode();
1516
1517   const gp_XY& P1=GetVertex(p1).Coord();
1518   const gp_XY& P2=GetVertex(p2).Coord();
1519   const gp_XY& P3=GetVertex(p3).Coord();
1520   gp_XY E1(P2); E1.Subtract(P1);
1521   gp_XY E2(P3); E2.Subtract(P2);
1522   gp_XY E3(P1); E3.Subtract(P3);
1523
1524   Standard_Real mode1=E1.SquareModulus();
1525   //Standard_Real dist=Sqrt(mode1);
1526   if (mode1<=EPSEPS) return Standard_False;
1527   Standard_Real v1=E1^(vert.Coord()-P1);
1528   Standard_Real distMin=(v1*v1)/mode1;
1529   edgOn=e1;
1530
1531   Standard_Real mode2=E2.SquareModulus();
1532   Standard_Real dist;
1533   //dist=Sqrt(mode2);
1534   if (mode2<=EPSEPS) return Standard_False;
1535   Standard_Real v2=E2^(vert.Coord()-P2);
1536   dist=(v2*v2)/mode2;
1537   if (dist<distMin) {
1538     edgOn=e2;
1539     distMin=dist;
1540   }
1541
1542   Standard_Real mode3=E3.SquareModulus();
1543   //dist=Sqrt(mode3);
1544   if (mode3<=EPSEPS) return Standard_False;
1545   Standard_Real v3=E3^(vert.Coord()-P3);
1546   dist=(v3*v3)/mode3;
1547   if (dist<distMin) {
1548     edgOn=e3;
1549     distMin=dist;
1550   }
1551
1552   if (distMin>EPSEPS) {
1553     Standard_Integer edf=edgOn;
1554     edgOn=0;
1555     if      (edf==e1 && edg1.Movability()!=MeshDS_Free) {
1556       if (v1<(mode1/5.)) edgOn=e1;
1557     }
1558     else if (edf==e2 && edg2.Movability()!=MeshDS_Free) {
1559       if (v2<(mode2/5.)) edgOn=e2;
1560     }
1561     else if (edf==e3 && edg3.Movability()!=MeshDS_Free) {
1562       if (v3<(mode3/5.)) edgOn=e3;
1563     }
1564   }
1565
1566   return (v1+v2+v3!=0. &&((v1>=0. && v2>=0. && v3>=0.) ||
1567                           (v1<=0. && v2<=0. && v3<=0.)));
1568 }
1569
1570 //=======================================================================
1571 //function : TriangleContaining
1572 //purpose  :
1573 //=======================================================================
1574 Standard_Integer MeshAlgo_Delaunay::TriangleContaining(const Vertex& vert)
1575 {
1576   const MeshDS_ListOfInteger& cirL=tCircles.Select(vert.Coord());
1577
1578   MeshDS_ListOfInteger::Iterator itT(cirL);
1579   Standard_Integer triPer=0;
1580   Standard_Integer edgon=0;
1581   for (; itT.More(); itT.Next()) {
1582     if (Contains(itT.Value(), vert, edgon)) {
1583       if (edgon==0) {
1584         triPer=itT.Value();
1585         break;
1586       }
1587       else if (GetEdge(edgon).Movability()==MeshDS_Free) {
1588         triPer=itT.Value();
1589         break;
1590       }
1591     }
1592   }
1593   return triPer;
1594 }