src/SMESHUtils/SMESH_FillHole.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 // File      : SMESH_FillHole.cxx
  23 // Created   : Tue Sep 26 15:11:17 2017
  24 // Author    : Edward AGAPOV (eap)
  25 //
  26
  27 #include "SMESH_MeshAlgos.hxx"
  28
  29 #include "SMESH_Comment.hxx"
  30
  31 #include "SMESH_TypeDefs.hxx"
  32 #include "SMDS_Mesh.hxx"
  33
  34 #include <Utils_SALOME_Exception.hxx>
  35
  36 #include <boost/intrusive/circular_list_algorithms.hpp>
  37 #include <boost/container/flat_map.hpp>
  38
  39 #include <Bnd_B3d.hxx>
  40
  41 namespace
  42 {
  43   bool isSmallAngle( double cos2 )
  44   {
  45     // cosine of min angle at which adjacent faces are considered overlapping
  46     const double theMinCos2 = 0.996 * 0.996; // ~5 degrees
  47     return ( cos2 > theMinCos2 );
  48   }
  49
  50   struct BEdge;
  51   typedef std::multimap< double, BEdge* >          TAngleMap;
  52   typedef std::map< const SMDS_MeshElement*, int > TFaceIndMap;
  53
  54   //--------------------------------------------------------------------------------
  55   /*!
  56    * \brief Edge of a free border
  57    */
  58   struct BEdge
  59   {
  60     const SMDS_MeshNode*    myNode1;
  61     const SMDS_MeshNode*    myNode2;
  62     const SMDS_MeshElement* myFace;   // face adjacent to the border
  63
  64     gp_XYZ                  myFaceNorm;
  65     gp_XYZ                  myDir;     // myNode1 -> myNode2
  66     double                  myDirCoef; // 1. or -1, to make myDir oriented as myNodes in myFace
  67     double                  myLength;  // between nodes
  68     double                  myAngleWithPrev; // between myDir and -myPrev->myDir
  69     double                  myMinMaxRatio; // of a possible triangle sides
  70     TAngleMap::iterator     myAngleMapPos;
  71     double                  myOverlapAngle;  // angle delta due to overlapping
  72     const SMDS_MeshNode*    myNode1Shift;    // nodes created to avoid overlapping of faces
  73     const SMDS_MeshNode*    myNode2Shift;
  74
  75     BEdge*                  myPrev; // neighbors in the border
  76     BEdge*                  myNext;
  77
  78     BEdge(): myNode1Shift(0), myNode2Shift(0) {}
  79     void   Init( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2,
  80                  const SMDS_MeshElement* f=0,
  81                  const SMDS_MeshNode* nf1=0, const SMDS_MeshNode* nf2=0 );
  82     void   ComputeAngle( bool reverseAngle = false );
  83     void   ShiftOverlapped( const SMDS_MeshNode*                  oppNode,
  84                             const TFaceIndMap&                    capFaceWithBordInd,
  85                             SMDS_Mesh&                            mesh,
  86                             std::vector<const SMDS_MeshElement*>& newFaces);
  87     void   MakeShiftfFaces( SMDS_Mesh&                            mesh,
  88                             std::vector<const SMDS_MeshElement*>& newFaces,
  89                             const bool                            isReverse );
  90     gp_XYZ GetInFaceDir() const { return myFaceNorm ^ myDir * myDirCoef; }
  91     double ShapeFactor()  const { return 0.5 * ( 1. - myMinMaxRatio ); }
  92     void   InsertSelf(TAngleMap& edgesByAngle, bool isReverseFaces, bool reBind, bool useOverlap )
  93     {
  94       if ( reBind ) edgesByAngle.erase( myAngleMapPos );
  95       double key = (( isReverseFaces ? 2 * M_PI - myAngleWithPrev : myAngleWithPrev )
  96                     + myOverlapAngle * useOverlap
  97                     + ShapeFactor() );
  98       myAngleMapPos = edgesByAngle.insert( std::make_pair( key, this ));
  99     }
 100
 101     // traits used by boost::intrusive::circular_list_algorithms
 102     typedef BEdge         node;
 103     typedef BEdge *       node_ptr;
 104     typedef const BEdge * const_node_ptr;
 105     static node_ptr get_next(const_node_ptr n)             {  return n->myNext;  }
 106     static void     set_next(node_ptr n, node_ptr next)    {  n->myNext = next;  }
 107     static node_ptr get_previous(const_node_ptr n)         {  return n->myPrev;  }
 108     static void     set_previous(node_ptr n, node_ptr prev){  n->myPrev = prev;  }
 109   };
 110
 111   //================================================================================
 112   /*!
 113    * \brief Initialize a border edge data
 114    */
 115   //================================================================================
 116
 117   void BEdge::Init( const SMDS_MeshNode*    n1,
 118                     const SMDS_MeshNode*    n2,
 119                     const SMDS_MeshElement* newFace, // new cap face
 120                     const SMDS_MeshNode*    nf1,
 121                     const SMDS_MeshNode*    nf2 )
 122   {
 123     myNode1  = n1;
 124     myNode2  = n2;
 125     myDir    = SMESH_NodeXYZ( n2 ) - SMESH_NodeXYZ( n1 );
 126     myLength = myDir.Modulus();
 127     if ( myLength > std::numeric_limits<double>::min() )
 128       myDir /= myLength;
 129
 130     myFace = newFace;
 131     if ( !myFace )
 132     {
 133       TIDSortedElemSet elemSet, avoidSet;
 134       int ind1, ind2;
 135       myFace = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet, &ind1, &ind2 );
 136       if ( !myFace )
 137         throw SALOME_Exception( SMESH_Comment("No face sharing nodes #")
 138                                 << myNode1->GetID() << " and #" << myNode2->GetID());
 139       avoidSet.insert( myFace );
 140       if ( SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet ))
 141         throw SALOME_Exception( SMESH_Comment("No free border between nodes #")
 142                                 << myNode1->GetID() << " and #" << myNode2->GetID());
 143
 144       myDirCoef = SMESH_MeshAlgos::IsRightOrder( myFace, myNode1, myNode2 ) ? 1. : -1.;
 145     }
 146
 147     if (! SMESH_MeshAlgos::FaceNormal( myFace, myFaceNorm, /*normalized=*/false ))
 148     {
 149       SMDS_ElemIteratorPtr fIt = myNode1->GetInverseElementIterator( SMDSAbs_Face );
 150       while ( fIt->more() )
 151         if ( SMESH_MeshAlgos::FaceNormal( fIt->next(), myFaceNorm, /*normalized=*/false ))
 152           break;
 153     }
 154
 155     if ( newFace )
 156     {
 157       myFace    = 0;
 158       myDirCoef = SMESH_MeshAlgos::IsRightOrder( newFace, nf1, nf2 ) ? 1. : -1.;
 159       if ( myPrev->myNode2 == n1 )
 160         myNode1Shift = myPrev->myNode2Shift;
 161       if ( myNext->myNode1 == n2 )
 162         myNode2Shift = myNext->myNode1Shift;
 163     }
 164     else if ( myDirCoef * myPrev->myDirCoef < 0 ) // different orientation of faces
 165     {
 166       myFaceNorm *= -1;
 167       myDirCoef  *= -1;
 168     }
 169   }
 170
 171   //================================================================================
 172   /*!
 173    * \brief Compute myAngleWithPrev
 174    */
 175   //================================================================================
 176
 177   void BEdge::ComputeAngle( bool theReverseAngle )
 178   {
 179     double dot = myDir.Dot( myPrev->myDir.Reversed() );
 180     if      ( dot >=  1 ) myAngleWithPrev = 0;
 181     else if ( dot <= -1 ) myAngleWithPrev = M_PI;
 182     else                  myAngleWithPrev = acos( dot );
 183
 184     bool isObtuse;
 185     gp_XYZ inFaceDirNew = myDir - myPrev->myDir;
 186     gp_XYZ   inFaceDir1 = myPrev->GetInFaceDir();
 187     gp_XYZ   inFaceDir2 = this->GetInFaceDir();
 188     double         dot1 = inFaceDirNew * inFaceDir1;
 189     double         dot2 = inFaceDirNew * inFaceDir2;
 190     bool     isOverlap1 = ( dot1 > 0 );
 191     bool     isOverlap2 = ( dot2 > 0 );
 192     if ( !myPrev->myFace )
 193       isObtuse = isOverlap1;
 194     else if  ( !myFace )
 195       isObtuse = isOverlap2;
 196     else
 197     {
 198       double dt1 = myDir.Dot( myPrev->myFaceNorm );
 199       double dt2 = myPrev->myDir.Dot( myFaceNorm );
 200       isObtuse = ( dt1 > 0 || dt2 < 0 ); // suppose face normals point outside the border
 201       if ( theReverseAngle )
 202         isObtuse = !isObtuse;
 203     }
 204     if ( isObtuse )
 205     {
 206       myAngleWithPrev = 2 * M_PI - myAngleWithPrev;
 207     }
 208
 209     // if ( ! isObtuse )
 210     //   isObtuse =
 211     //     isSmallAngle( 1 - myDir.CrossSquareMagnitude( myPrev->myDir )); // edges co-directed
 212
 213     myOverlapAngle = 0;
 214     //if ( !isObtuse )
 215     {
 216       // check if myFace and a triangle built on this and prev edges overlap
 217       if ( isOverlap1 )
 218       {
 219         double cos2 = dot1 * dot1 / inFaceDirNew.SquareModulus() / inFaceDir1.SquareModulus();
 220         myOverlapAngle += 1. * M_PI * cos2;
 221       }
 222       if ( isOverlap2 )
 223       {
 224         double cos2 = dot2 * dot2 / inFaceDirNew.SquareModulus() / inFaceDir2.SquareModulus();
 225         myOverlapAngle += 1. * M_PI * cos2;
 226       }
 227     }
 228
 229     {
 230       double len3 = SMESH_NodeXYZ( myPrev->myNode1 ).Distance( myNode2 );
 231       double minLen = Min( myLength, Min( myPrev->myLength, len3 ));
 232       double maxLen = Max( myLength, Max( myPrev->myLength, len3 ));
 233       myMinMaxRatio = minLen / maxLen;
 234     }
 235   }
 236
 237   //================================================================================
 238   /*!
 239    * \brief Check if myFace is overlapped by a triangle formed by myNode's and a
 240    *        given node. If so, create shifted nodes to avoid overlapping
 241    */
 242   //================================================================================
 243
 244   void BEdge::ShiftOverlapped( const SMDS_MeshNode*                  theOppNode,
 245                                const TFaceIndMap&                    theCapFaceWithBordInd,
 246                                SMDS_Mesh&                            theMesh,
 247                                std::vector<const SMDS_MeshElement*>& theNewFaces )
 248   {
 249     if ( myNode1Shift && myNode2Shift )
 250       return;
 251
 252     gp_XYZ inNewFaceDir = SMESH_NodeXYZ( theOppNode ) - SMESH_NodeXYZ( myNode1 );
 253     double          dot = inNewFaceDir.Dot( myFaceNorm );
 254     double         cos2 = dot * dot / myFaceNorm.SquareModulus() / inNewFaceDir.SquareModulus();
 255     bool      isOverlap = ( isSmallAngle( 1 - cos2 ) && GetInFaceDir() * inNewFaceDir > 0 );
 256
 257     if ( isOverlap )
 258     {
 259       gp_XYZ shift = myFaceNorm / myLength / 4;
 260       if ( myFace )
 261         shift.Reverse();
 262       if ( !myNode1Shift )
 263       {
 264         gp_XYZ p = SMESH_NodeXYZ( myNode1 ) + shift;
 265         myNode1Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
 266         myPrev->myNode2Shift = myNode1Shift;
 267       }
 268       if ( !myNode2Shift )
 269       {
 270         gp_XYZ p = SMESH_NodeXYZ( myNode2 ) + shift;
 271         myNode2Shift = theMesh.AddNode( p.X(), p.Y(), p.Z() );
 272         myNext->myNode1Shift = myNode2Shift;
 273       }
 274
 275       // MakeShiftfFaces() for already created cap faces
 276       for ( int is2nd = 0; is2nd < 2; ++is2nd )
 277       {
 278         const SMDS_MeshNode* ns = is2nd ? myNode2Shift : myNode1Shift;
 279         const SMDS_MeshNode* n  = is2nd ? myNode2 : myNode1;
 280         if ( !ns ) continue;
 281
 282         SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
 283         while ( fIt->more() )
 284         {
 285           const SMDS_MeshElement* f = fIt->next();
 286           if ( !f->isMarked() ) continue;
 287
 288           TFaceIndMap::const_iterator f2i = theCapFaceWithBordInd.find( f );
 289           if ( f2i == theCapFaceWithBordInd.end() )
 290             continue;
 291           const SMDS_MeshNode* nf1 = f->GetNode(  f2i->second );
 292           const SMDS_MeshNode* nf2 = f->GetNode(( f2i->second+1 ) % f->NbNodes() );
 293           if ( nf1 == n || nf2 == n )
 294           {
 295             BEdge tmpE;
 296             tmpE.myPrev = tmpE.myNext = this;
 297             tmpE.Init( nf1, nf2, f, nf1, nf2 );
 298             if ( !tmpE.myNode1Shift && !tmpE.myNode2Shift )
 299               tmpE.Init( nf2, nf1, f, nf2, nf1 );
 300             tmpE.myFace = f;
 301             tmpE.MakeShiftfFaces( theMesh, theNewFaces, tmpE.myDirCoef < 0 );
 302           }
 303           std::vector< const SMDS_MeshNode* > nodes( f->begin_nodes(), f->end_nodes() );
 304           nodes[ f->GetNodeIndex( n ) ] = ns;
 305           theMesh.ChangeElementNodes( f, &nodes[0], nodes.size() );
 306         }
 307       }
 308     }
 309   }
 310
 311   //================================================================================
 312   /*!
 313    * \brief Create a triangle
 314    */
 315   //================================================================================
 316
 317   const SMDS_MeshElement* MakeTria( SMDS_Mesh&           mesh,
 318                                     const SMDS_MeshNode* n1,
 319                                     const SMDS_MeshNode* n2,
 320                                     const SMDS_MeshNode* n3,
 321                                     const bool           isReverse )
 322   {
 323     if ( isReverse )
 324       return mesh.AddFace( n1, n3, n2 );
 325     return mesh.AddFace( n1, n2, n3 );
 326   }
 327
 328   //================================================================================
 329   /*!
 330    * \brief Create a quadrangle
 331    */
 332   //================================================================================
 333
 334   // const SMDS_MeshElement* MakeQuad( SMDS_Mesh&           mesh,
 335   //                                   const SMDS_MeshNode* n1,
 336   //                                   const SMDS_MeshNode* n2,
 337   //                                   const SMDS_MeshNode* n3,
 338   //                                   const SMDS_MeshNode* n4,
 339   //                                   const bool           isReverse )
 340   // {
 341   //   if ( isReverse )
 342   //     return mesh.AddFace( n4, n3, n2, n1 );
 343   //   return mesh.AddFace( n1, n2, n3, n4 );
 344   // }
 345
 346   //================================================================================
 347   /*!
 348    * \brief Create faces on myNode* and myNode*Shift
 349    */
 350   //================================================================================
 351
 352   void BEdge::MakeShiftfFaces(SMDS_Mesh&                            mesh,
 353                               std::vector<const SMDS_MeshElement*>& newFaces,
 354                               const bool                            isReverse )
 355   {
 356     if ( !myFace )
 357       return;
 358     if ( myNode1Shift && myNode2Shift )
 359     {
 360       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
 361       newFaces.push_back( MakeTria( mesh, myNode1, myNode2Shift, myNode1Shift, isReverse ));
 362     }
 363     else if ( myNode1Shift )
 364     {
 365       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode1Shift, isReverse ));
 366     }
 367     else if ( myNode2Shift )
 368     {
 369       newFaces.push_back( MakeTria( mesh, myNode1, myNode2, myNode2Shift, isReverse ));
 370     }
 371   }
 372
 373 } // namespace
 374
 375 //================================================================================
 376 /*!
 377  * \brief Fill with 2D elements a hole defined by a TFreeBorder
 378  */
 379 //================================================================================
 380
 381 void SMESH_MeshAlgos::FillHole(const SMESH_MeshAlgos::TFreeBorder &  theFreeBorder,
 382                                SMDS_Mesh&                            theMesh,
 383                                std::vector<const SMDS_MeshElement*>& theNewFaces)
 384 {
 385   if ( theFreeBorder.size() < 4 ||                // at least 3 nodes
 386        theFreeBorder[0] != theFreeBorder.back() ) // the hole must be closed
 387     return;
 388
 389   // prepare data of the border
 390
 391   ObjectPool< BEdge > edgeAllocator;
 392   boost::intrusive::circular_list_algorithms< BEdge > circularList;
 393   BEdge* edge;
 394   BEdge* edge0 = edgeAllocator.getNew();
 395   BEdge* edgePrev = edge0;
 396   circularList.init_header( edge0 );
 397   edge0->Init( theFreeBorder[0], theFreeBorder[1], 0 );
 398   Bnd_B3d box;
 399   box.Add( SMESH_NodeXYZ( edge0->myNode1 ));
 400   for ( size_t i = 2; i < theFreeBorder.size(); ++i )
 401   {
 402     edge = edgeAllocator.getNew();
 403     circularList.link_after( edgePrev, edge );
 404     edge->Init( theFreeBorder[i-1], theFreeBorder[i] );
 405     edge->ComputeAngle();
 406     edgePrev = edge;
 407     box.Add( SMESH_NodeXYZ( edge->myNode1 ));
 408   }
 409   edge0->ComputeAngle();
 410
 411   // check if face normals point outside the border
 412
 413   gp_XYZ hSize = 0.5 * ( box.CornerMax() - box.CornerMin() );
 414   const double hDelta = 1e-6 * hSize.Modulus();
 415   hSize -= gp_XYZ( hDelta, hDelta, hDelta );
 416   if ( hSize.X() < 0 ) hSize.SetX(hDelta);
 417   if ( hSize.Y() < 0 ) hSize.SetY(hDelta);
 418   if ( hSize.Z() < 0 ) hSize.SetZ(hDelta);
 419   box.SetHSize( hSize ); // decrease the box by hDelta
 420
 421   size_t nbEdges = theFreeBorder.size() - 1;
 422   edge = edge0;
 423   int nbRev = 0, nbFrw = 0;
 424   double angTol = M_PI - ( nbEdges - 2 ) * M_PI / nbEdges, sumDirCoeff = 0;
 425   for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
 426   {
 427     if ( box.IsOut( SMESH_NodeXYZ( edge->myNode1 )) &&
 428          edge->myOverlapAngle < 0.1 * M_PI )
 429     {
 430       nbRev += edge->myAngleWithPrev > M_PI + angTol;
 431       nbFrw += edge->myAngleWithPrev < M_PI - angTol;
 432     }
 433     sumDirCoeff += edge->myDirCoef;
 434
 435     // unmark all adjacent faces, new faces will be marked
 436     SMDS_ElemIteratorPtr fIt = edge->myNode1->GetInverseElementIterator( SMDSAbs_Face );
 437     while ( fIt->more() )
 438       fIt->next()->setIsMarked( false );
 439   }
 440   bool isReverseAngle = ( nbRev > nbFrw ); // true == face normals point inside the border
 441   //std::cout << "nbRev="<< nbRev << ", nbFrw="<< nbFrw<<std::endl;
 442
 443   // sort border edges by myAngleWithPrev
 444
 445   TAngleMap edgesByAngle;
 446   bool useOverlap = true; // to add BEdge.myOverlapAngle when filling edgesByAngle
 447   edge = edge0;
 448   for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
 449     edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
 450
 451   // create triangles to fill the hole
 452
 453   //compare order of nodes in the edges with their order in faces
 454   bool isReverse = sumDirCoeff > 0.5 * nbEdges;
 455
 456   // faces filling the hole (cap faces) and indices of border edges in them
 457   TFaceIndMap capFaceWithBordInd;
 458
 459   theNewFaces.reserve( nbEdges - 2 );
 460   std::vector< const SMDS_MeshNode* > nodes(3);
 461   while ( edgesByAngle.size() > 2 )
 462   {
 463     TAngleMap::iterator a2e = edgesByAngle.begin();
 464     edge = a2e->second;
 465     if ( useOverlap &&
 466          a2e->first - edge->ShapeFactor() > M_PI - angTol ) // all new triangles need shift
 467     {
 468       // re-sort the edges w/o overlap consideration
 469       useOverlap = false;
 470       nbEdges = edgesByAngle.size();
 471       edgesByAngle.clear();
 472       for ( size_t i = 0; i < nbEdges; ++i, edge = edge->myNext )
 473         edge->InsertSelf( edgesByAngle, isReverseAngle, /*reBind=*/false, useOverlap );
 474       a2e = edgesByAngle.begin();
 475     }
 476     edge     = a2e->second;
 477     edgePrev = edge->myPrev;
 478
 479     // create shift nodes and faces
 480     edgePrev->ShiftOverlapped( edge->myNode2, capFaceWithBordInd, theMesh, theNewFaces );
 481     edge->ShiftOverlapped( edgePrev->myNode1, capFaceWithBordInd, theMesh, theNewFaces );
 482     edge    ->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
 483     edgePrev->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
 484
 485     // make a cap face
 486     //nodes.resize( 3 );
 487     nodes[0] = edgePrev->myNode1Shift ? edgePrev->myNode1Shift : edgePrev->myNode1;
 488     nodes[1] = edgePrev->myNode2Shift ? edgePrev->myNode2Shift : edgePrev->myNode2;
 489     nodes[2] = edge->myNode2Shift     ? edge->myNode2Shift     : edge->myNode2;
 490     theNewFaces.push_back( MakeTria( theMesh, nodes[0], nodes[1], nodes[2], isReverse ));
 491     // std::cout << nodes[1]->GetID() << " "  << nodes[0]->GetID() << " "  << nodes[2]->GetID()
 492     //           << " " << edge->myAngleWithPrev << std::endl;
 493
 494     // remember a border edge within the new cap face
 495     theNewFaces.back()->setIsMarked( true );
 496     if ( edgePrev->myFace )
 497       capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), isReverse ? 2 : 0 ));
 498     if ( edge->myFace )
 499       capFaceWithBordInd.insert( std::make_pair( theNewFaces.back(), 1 ));
 500
 501     // remove edgePrev from the list and update <edge>
 502     edgesByAngle.erase( edgePrev->myAngleMapPos );
 503     circularList.unlink( edgePrev ); // remove edgePrev from the border
 504
 505     edge->Init( edgePrev->myNode1, edge->myNode2, theNewFaces.back(), nodes[0], nodes[2] );
 506     edge->ComputeAngle( isReverseAngle );
 507     edge->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
 508     edge->myNext->ComputeAngle( isReverseAngle );
 509     edge->myNext->InsertSelf( edgesByAngle, /*isReverse=*/false, /*reBind=*/true, useOverlap );
 510     // std::cout << "A " << edge->myNode1->GetID() << " " << edge->myAngleWithPrev
 511     //           << " " << edge->myNext->myNode1->GetID() << " " << edge->myNext->myAngleWithPrev
 512     //           << std::endl;
 513   }
 514   edge = edgesByAngle.begin()->second;
 515   edge->        MakeShiftfFaces( theMesh, theNewFaces, isReverse );
 516   edge->myNext->MakeShiftfFaces( theMesh, theNewFaces, isReverse );
 517 }