src/BLSURFPlugin/BLSURFPlugin_EnforcedMesh1D.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // File      : BLSURFPlugin_EnforcedMesh1D.cxx
  20 // Author    : Edward AGAPOV (OCC)
  21
  22 #include "BLSURFPlugin_EnforcedMesh1D.hxx"
  23
  24 #include "BLSURFPlugin_BLSURF.hxx"
  25
  26 #include <SMDS_IteratorOnIterators.hxx>
  27 #include <SMDS_MeshEdge.hxx>
  28 #include <SMDS_MeshGroup.hxx>
  29 #include <SMESHDS_Group.hxx>
  30 #include <SMESHDS_Mesh.hxx>
  31 #include <SMESH_Group.hxx>
  32 #include <SMESH_Mesh.hxx>
  33 #include <SMESH_MeshAlgos.hxx>
  34 #include <SMESH_MeshEditor.hxx>
  35 #include <SMESH_MesherHelper.hxx>
  36 #include <SMESH_TypeDefs.hxx>
  37
  38 #include <BRepAdaptor_Curve.hxx>
  39 #include <BRepBuilderAPI_MakeVertex.hxx>
  40 #include <BRep_Tool.hxx>
  41 #include <Extrema_ExtCC.hxx>
  42 #include <Extrema_ExtElC.hxx>
  43 #include <Geom_Line.hxx>
  44 #include <ShapeAnalysis_Curve.hxx>
  45 #include <Geom2d_Line.hxx>
  46
  47 // allow range iteration on NCollection_IndexedMap
  48 template < class IMAP >
  49 typename IMAP::const_iterator begin( const IMAP &  m ) { return m.cbegin(); }
  50 template < class IMAP >
  51 typename IMAP::const_iterator end( const IMAP &  m ) { return m.cend(); }
  52
  53
  54 namespace
  55 {
  56   //================================================================================
  57   /*!
  58    * \brief Look for a node coincident with a given nodes among end nodes of a segment
  59    *        and among its internal nodes
  60    *  \param [in] p - the epoint
  61    *  \param [in] tol - 3D tolarace
  62    *  \param [in] segment - the segment
  63    *  \param [in] segInternalNodes - map of segment internal nodes
  64    *  \param [out] index - return index of found internal node; -1 if an end node is found
  65    *  \return const SMDS_MeshNode* - found node
  66    */
  67   //================================================================================
  68
  69   const SMDS_MeshNode* findNode( const gp_Pnt&                              p,
  70                                  const double                               tol,
  71                                  const SMDS_MeshElement*                    segment,
  72                                  BLSURFPlugin_EnforcedMesh1D::TNodesOnSeg & segInternalNodes,
  73                                  int &                                      index )
  74   {
  75
  76     SMESH_NodeXYZ node0 = segment->GetNode(0);
  77     if ( p.IsEqual( node0, tol ))
  78       return index = -1, node0.Node();
  79     SMESH_NodeXYZ node1 = segment->GetNode(1);
  80     if ( p.IsEqual( node1, tol ))
  81       return index = -1, node1.Node();
  82
  83     auto seg2nodes = segInternalNodes.find( segment );
  84     if ( seg2nodes != segInternalNodes.end() )
  85     {
  86       const std::vector< const SMDS_MeshNode* >& nodes = seg2nodes->second;
  87       for ( size_t i = 0; i < nodes.size(); ++i )
  88         if ( p.IsEqual( SMESH_NodeXYZ( nodes[i] ), tol ))
  89         {
  90           index = (int) i;
  91           return nodes[i];
  92         }
  93     }
  94     return nullptr;
  95   }
  96
  97   //================================================================================
  98   /*!
  99    * \brief Orient segments to correspond to order of nodes in a branch
 100    *  \param [in] braSegs - segments of the branch
 101    *  \param [in] braNodes - nodes of the branch
 102    *  \param [in] nodeIndex - index of a node of the branch
 103    *  \param [inout] mesh - mesh holding the nodes and segments
 104    *
 105    *
 106    */
 107   //================================================================================
 108
 109   void orientSegments( const std::vector< const SMDS_MeshElement* > & braSegs,
 110                        const std::vector< const SMDS_MeshNode* > &    braNodes,
 111                        const size_t                                   nodeIndex,
 112                        SMESH_Mesh*                                    mesh )
 113   {
 114     const SMDS_MeshElement* toReverse[2] = { nullptr, nullptr };
 115
 116     if ( nodeIndex > 0 &&
 117          braSegs[ nodeIndex - 1 ]->GetNode(1) != braNodes[ nodeIndex ])
 118       toReverse[ 0 ] = braSegs[ nodeIndex - 1 ];
 119
 120     if ( nodeIndex < braSegs.size() &&
 121          braSegs[ nodeIndex ]->GetNode(0) != braNodes[ nodeIndex ])
 122       toReverse[ bool( toReverse[0]) ] = braSegs[ nodeIndex ];
 123
 124     if ( !toReverse[0] )
 125       return;
 126
 127     SMESH_MeshEditor editor( mesh );
 128     for ( int i = 0; i < 2; ++i )
 129       if ( toReverse[ i ])
 130         editor.Reorient( toReverse[ i ]);
 131   }
 132
 133 } // namespace
 134
 135 //================================================================================
 136 /*!
 137  * \brief Create enforced mesh segments in a mesh
 138  *  \param [in] helper - contains the mesh and the shape to compute
 139  *  \param [in] hyp - hypothesis containing data of enforced meshes
 140  */
 141 //================================================================================
 142
 143 BLSURFPlugin_EnforcedMesh1D::BLSURFPlugin_EnforcedMesh1D( SMESH_MesherHelper&            helper,
 144                                                           const BLSURFPlugin_Hypothesis* hyp )
 145   : _mesh ( helper.GetMesh() ),
 146     _shape( helper.GetSubShape() ),
 147     _helper( *_mesh ),
 148     _isQuadratic( helper.GetIsQuadratic() )
 149 {
 150   if ( !hyp || !_mesh || hyp->GetEnforcedMeshes().empty() )
 151     return;
 152
 153   _tol      = 2 * BRep_Tool::MaxTolerance( _shape, TopAbs_VERTEX );
 154   _segTag   = 1000 + helper.Count( _shape, TopAbs_EDGE, /*ignoreSame=*/false );
 155   _segTag0  = _segTag;
 156   _nodeTag0 = 1000 + helper.Count( _shape, TopAbs_VERTEX, /*ignoreSame=*/false );
 157
 158   for ( const BLSURFPlugin_Hypothesis::EnforcedMesh& enfMesh : hyp->GetEnforcedMeshes() )
 159   {
 160     copyEnforcedMesh( enfMesh, hyp, _shape );
 161   }
 162
 163   for ( const SMDS_MeshElement* segment : _segmentsOnSeveralShapes )
 164   {
 165     splitSegmentOnSeveralShapes( segment );
 166   }
 167
 168   for ( const TopoDS_Shape & face : _faces )
 169   {
 170     splitSelfIntersectingSegments( face );
 171   }
 172
 173   for ( TEdge2Nodes::Iterator e2nn( _nodesOnEdge ); e2nn.More(); e2nn.Next() )
 174   {
 175     splitEdgeByNodes( e2nn.Key(), e2nn.Value() );
 176   }
 177 }
 178
 179 //================================================================================
 180 /*!
 181  * \brief Convert enforced segments to quadratic
 182  */
 183 //================================================================================
 184
 185 BLSURFPlugin_EnforcedMesh1D::~BLSURFPlugin_EnforcedMesh1D()
 186 {
 187   if ( !_isQuadratic )
 188     return;
 189
 190   _helper.SetIsQuadratic( true );
 191
 192   SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
 193   std::vector<const SMDS_MeshNode *>    nodes(2);
 194   std::vector<const SMDS_MeshElement *> trias;
 195
 196   SMDS_MeshGroup* group = nullptr;
 197   if ( !_name2Group.empty() )
 198     group = _name2Group.begin()->second;
 199
 200   for ( const TopoDS_Shape& f : _faces )
 201     if ( HasSegmentsOnFace( TopoDS::Face( f )))
 202       while ( _segIterator->more() )
 203       {
 204         const SMDS_MeshElement* segment = _segIterator->next();
 205         if ( segment->GetType() != SMDSAbs_Edge )
 206           break;
 207
 208         nodes[0] = segment->GetNode(0);
 209         nodes[1] = segment->GetNode(1);
 210         if ( meshDS->GetElementsByNodes( nodes, trias, SMDSAbs_Face ))
 211         {
 212           if ( !group || !group->Contains( segment ))
 213           {
 214             SMDS_MeshGroup* otherGroup = nullptr;
 215             if ( _name2Group.size() > bool( group ))
 216               for ( auto & n2g : _name2Group )
 217                 if ( n2g.second != group && n2g.second->Contains( segment ))
 218                 {
 219                   otherGroup = n2g.second;
 220                   break;
 221                 }
 222             group = otherGroup;
 223           }
 224
 225           _helper.AddTLinks( static_cast< const SMDS_MeshFace* >( trias[0] ));
 226
 227           meshDS->RemoveFreeElement( segment, meshDS->MeshElements( f ), /*fromGroup=*/false );
 228
 229           SMDS_MeshElement* quadSegment = _helper.AddEdge( nodes[0], nodes[1],
 230                                                            /*id=*/0, /*force3d=*/false );
 231           if ( group && segment != quadSegment )
 232           {
 233             group->Remove( segment );
 234             group->Add( quadSegment );
 235           }
 236         }
 237       }
 238   return;
 239 }
 240
 241 //================================================================================
 242 /*!
 243  * \brief Return EDGEs resulted from division of FACE boundary by enforced segments
 244  *  \param [in] edge - possibly divided EDGE
 245  *  \param [out] splits - split EDGEs
 246  *  \return bool - true if the EDGE is split
 247  */
 248 //================================================================================
 249
 250 bool BLSURFPlugin_EnforcedMesh1D::GetSplitsOfEdge( const TopoDS_Edge&           edge,
 251                                                    std::vector< TopoDS_Edge > & splits,
 252                                                    TopTools_IndexedMapOfShape & edgeTags )
 253 {
 254   std::vector< TopoDS_Edge > * splitsInMap = _edgeSplitsOfEdge.ChangeSeek( edge );
 255   if ( !splitsInMap )
 256     return false;
 257
 258   splits.insert( splits.end(), splitsInMap->begin(), splitsInMap->end() );
 259
 260   int eTag = edgeTags.Add( edge );
 261
 262   size_t index = splits.size() - 1;
 263   for ( size_t i = 0; i < splitsInMap->size(); ++i, --index )
 264   {
 265     int splitTag = edgeTags.Add( splits[ index ]);
 266     _splitTags2EdgeTag[ splitTag ] = eTag;
 267
 268     if ( edge.Orientation() == TopAbs_REVERSED )
 269       splits[ index ].Reverse();
 270   }
 271
 272   return true;
 273 }
 274
 275
 276 //================================================================================
 277 /*!
 278  * \brief Add 1D elements to the mesh
 279  */
 280 //================================================================================
 281
 282 void BLSURFPlugin_EnforcedMesh1D::
 283
 284 copyEnforcedMesh( const BLSURFPlugin_Hypothesis::EnforcedMesh& theEnfMesh,
 285                   const BLSURFPlugin_Hypothesis*               theHyp,
 286                   const TopoDS_Shape&                          theShape)
 287 {
 288   SMESH_Mesh* mesh1D;
 289   SMDS_ElemIteratorPtr segIt = theHyp->GetEnforcedSegments( theEnfMesh, mesh1D );
 290   if ( !segIt->more() )
 291     return;
 292
 293   // setup predicates to detect nodes on FACE boundary
 294   setupPredicates( theShape );
 295
 296   SMDS_MeshGroup* group = getGroup( theEnfMesh._groupName );
 297   SMESHDS_Mesh*  meshDS = _helper.GetMeshDS();
 298
 299   // get ordered nodes and segments of theEnfMesh
 300   SMESH_MeshAlgos::TElemGroupVector edgeBranches;
 301   SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
 302   SMESH_MeshAlgos::Get1DBranches( segIt, edgeBranches, nodeBranches );
 303
 304
 305   // Copy nodes and segments from an enforced mesh to my mesh
 306
 307   TopoDS_Shape vertex, edge;
 308
 309   // first treat ends of branches that can be shared by branches
 310   for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
 311   {
 312     std::vector< const SMDS_MeshNode* > &   braNodes = nodeBranches[ iB ];
 313     std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
 314
 315     for ( int isLast = 0; isLast < 2; ++isLast )
 316     {
 317       const SMDS_MeshNode* enfNode = isLast ? braNodes.back() : braNodes[0];
 318       if ( meshDS->Contains( enfNode ))
 319         continue; // already in my mesh
 320
 321       const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
 322       if ( !newNode )
 323         orientSegments( braSegs, braNodes, isLast ? 0 : braNodes.size() - 1, mesh1D );
 324
 325       // replace enfNode at branch ends by newNode
 326       SMESH_NodeXYZ enfPnt( newNode ? newNode : enfNode );
 327       for ( std::vector< const SMDS_MeshNode* > & braNodes : nodeBranches )
 328       {
 329         for ( int isLast = 0; isLast < 2; ++isLast )
 330         {
 331           const SMDS_MeshNode* & endNode = isLast ? braNodes.back() : braNodes[0];
 332           if ( endNode == enfNode || enfPnt.SquareDistance( endNode ) < _tol*_tol )
 333             endNode = newNode;
 334         }
 335       }
 336       continue;
 337     }  // loop on branch ends
 338   } // loop on nodeBranches
 339
 340   // copy nodes and segments
 341
 342   for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
 343   {
 344     std::vector< const SMDS_MeshNode* > &   braNodes = nodeBranches[ iB ];
 345     std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
 346
 347     // copy nodes of the branch
 348     for ( size_t i = 0; i < braNodes.size(); ++i )
 349     {
 350       const SMDS_MeshNode* & enfNode = braNodes[ i ];
 351       const SMDS_MeshNode*   newNode = copyEnforcedNode( enfNode );
 352
 353       if ( !newNode ) // orient segments to be able to get enforced not projected node
 354         orientSegments( braSegs, braNodes, i, mesh1D );
 355
 356       enfNode = newNode;
 357     }
 358
 359     // copy segments of the branch
 360     for ( size_t i = 0; i < braSegs.size(); ++i )
 361     {
 362       //braSegs[ i ] = nullptr;
 363
 364       const SMDS_MeshNode* node0 = braNodes[ i     ];
 365       const SMDS_MeshNode* node1 = braNodes[ i + 1 ];
 366       if ( !node0 && !node1 )
 367         continue;
 368
 369       TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
 370       TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
 371       if ( shape0.IsNull() && shape1.IsNull() )
 372         continue;
 373
 374       if ( !node0 && shape1.ShapeType() != TopAbs_FACE )
 375         continue;
 376       if ( !node1 && shape0.ShapeType() != TopAbs_FACE )
 377         continue;
 378
 379       if ( !node0 || !node1 ) // create missing node at location of enforced node projected nowhere
 380       {
 381         SMESH_NodeXYZ  pn = braSegs[i]->GetNode( !node1 );
 382         ( node0 ? node1 : node0 ) = _helper.AddNode( pn->X(), pn->Y(), pn->Z() );
 383       }
 384
 385       SMDS_MeshEdge* newSeg = _helper.AddEdge( node0, node1 );
 386       if ( group )
 387         group->Add( newSeg );
 388       braSegs[ i ] = newSeg;
 389
 390       // check if the both nodes are on the same FACE
 391       TopoDS_Shape face = shape0;
 392       if ( !shape0.IsSame( shape1 ) && !shape0.IsNull() && !shape1.IsNull() )
 393       {
 394         if ( shape0.ShapeType() == TopAbs_FACE &&
 395              _helper.IsSubShape( shape1, shape0 ))
 396         {
 397           face = shape0;
 398         }
 399         else if ( shape1.ShapeType() == TopAbs_FACE &&
 400                   _helper.IsSubShape( shape0, shape1 ))
 401         {
 402           face = shape1;
 403         }
 404         else // try to find a FACE by projecting a segment middle point
 405         {
 406           face.Nullify();
 407           gp_Pnt middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
 408           //BLSURFPlugin_BLSURF::projectionPoint projPnt =
 409           BLSURFPlugin_BLSURF::getProjectionPoint( TopoDS::Face( face ), middlePnt );
 410
 411           if ( !face.IsNull() &&
 412                ( !_helper.IsSubShape( shape0, face ) ||
 413                  !_helper.IsSubShape( shape1, face ) ))
 414             face.Nullify();
 415         }
 416       }
 417       if ( !face.IsNull() && face.ShapeType() == TopAbs_FACE )
 418       {
 419         meshDS->SetMeshElementOnShape( newSeg, face );
 420         _faces.Add( face );
 421       }
 422
 423       if ( face.IsNull() || shape0.IsNull() || shape1.IsNull() )
 424       {
 425         _segmentsOnSeveralShapes.push_back( newSeg );
 426       }
 427
 428     } // loop on branch segments
 429     continue;
 430   } // loop on branches
 431
 432   return;
 433 }
 434
 435 //================================================================================
 436 /*!
 437  * \brief Create a copy of a node of enforced mesh in my mesh
 438  *  \param [in] enfNode - enforced node
 439  *  \return const SMDS_MeshNode* - a node in my mesh
 440  */
 441 //================================================================================
 442
 443 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::copyEnforcedNode( const SMDS_MeshNode* enfNode )
 444 {
 445   SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
 446
 447   if ( !enfNode || meshDS->Contains( enfNode ))
 448     return enfNode; // already in my mesh
 449
 450   SMESH_NodeXYZ enfPnt = enfNode;
 451
 452   const SMDS_MeshNode* newNode = nullptr;
 453
 454   // check if enfNode is on VERTEX
 455   TopoDS_Vertex vertex;
 456   if ( _onVertexPredicate->IsSatisfy( enfNode, &vertex ))
 457   {
 458     _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 459     newNode = SMESH_Algo::VertexNode( vertex, meshDS );
 460   }
 461
 462   // check if enfNode is on EDGE
 463   bool setNodeOnEdge = false;
 464   TopoDS_Edge edge;
 465   if ( !newNode )
 466   {
 467     setNodeOnEdge = _onEdgePredicate->IsSatisfy( enfNode, &edge );
 468     if ( setNodeOnEdge )
 469     {
 470       newNode = findNodeOnEdge( enfPnt, edge );
 471       setNodeOnEdge = !newNode;
 472     }
 473   }
 474
 475   // create a new node and set it on FACE
 476   if ( !newNode )
 477   {
 478     TopoDS_Face face;
 479     BLSURFPlugin_BLSURF::projectionPoint projPnt =
 480       BLSURFPlugin_BLSURF::getProjectionPoint( face, enfPnt, /*allowStateON=*/true );
 481     if ( face.IsNull() ) return newNode;
 482
 483     if ( projPnt.state == TopAbs_ON )
 484     {
 485       SMDS_MeshNode* projNode = const_cast< SMDS_MeshNode* >( enfNode );
 486       projNode->setXYZ( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
 487       vertex.Nullify();
 488       edge.Nullify();
 489       if ( _onVertexPredicate->IsSatisfy( projNode, &vertex ))
 490       {
 491         _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 492         newNode = SMESH_Algo::VertexNode( vertex, meshDS );
 493       }
 494       else if (( setNodeOnEdge = _onEdgePredicate->IsSatisfy( projNode, &edge )))
 495       {
 496         newNode = findNodeOnEdge( projPnt.xyz, edge );
 497         setNodeOnEdge = !newNode;
 498       }
 499       projNode->setXYZ( enfPnt.X(), enfPnt.Y(), enfPnt.Z() );
 500     }
 501
 502     if ( !newNode )
 503       newNode = meshDS->AddNode( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
 504
 505     if ( vertex.IsNull() && edge.IsNull() )
 506       meshDS->SetNodeOnFace( newNode, face, projPnt.uv.X(), projPnt.uv.Y() );
 507   }
 508
 509   // set the new node on EDGE
 510   if ( newNode && setNodeOnEdge )
 511   {
 512     addNodeOnEdge( newNode, edge );
 513   }
 514
 515   return newNode;
 516 }
 517
 518 //================================================================================
 519 /*!
 520  * \brief Split a segment whose nodes are on different FACEs into smaller parts
 521  *        lying each on one FACE
 522  */
 523 //================================================================================
 524
 525 void BLSURFPlugin_EnforcedMesh1D::splitSegmentOnSeveralShapes( const SMDS_MeshElement* segment )
 526 {
 527   const SMDS_MeshNode* node0 = segment->GetNode(0);
 528   const SMDS_MeshNode* node1 = segment->GetNode(1);
 529   SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
 530
 531   TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
 532   TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
 533
 534   if ( shape0.IsNull() )
 535   {
 536     std::swap( node0, node1 );
 537     shape0 = shape1;
 538     shape1.Nullify();
 539   }
 540
 541   gp_XYZ   xyz0 = SMESH_NodeXYZ( node0 );
 542   gp_XYZ   xyz1 = SMESH_NodeXYZ( node1 );
 543   gp_XYZ segDir = ( xyz1 - xyz0 ).Normalized();
 544
 545   //std::map< double, const SMDS_MeshNode* > mediumNodes; // nodes splitting the segment
 546
 547   while ( !shape0.IsSame( shape1 )) // move along the segment till shape1
 548   {
 549     if ( shape0.ShapeType() == TopAbs_FACE ) // make a node on an EDGE of the FACE
 550     {                                        // ----------------------------------
 551       TopoDS_Edge edge;
 552       double paramOnE;
 553       gp_Pnt edgeIntPnt = getEdgeIntersection( shape0, xyz0, xyz1, edge, paramOnE );
 554       if ( edge.IsNull() )
 555         break;
 556
 557       // check if edgeIntPnt in on VERTEX
 558       TopoDS_Vertex vertex;
 559       for ( int iV = 0; iV < 2 &&  vertex.IsNull(); ++iV )
 560       {
 561         TopoDS_Vertex v = _helper.IthVertex( iV, edge );
 562         if ( edgeIntPnt.SquareDistance( BRep_Tool::Pnt( v )) < _tol *_tol )
 563           vertex = v;
 564       }
 565
 566       // make a node on the EDGE
 567       const SMDS_MeshNode* nodeOnEdge = nullptr;
 568       if ( vertex.IsNull() )
 569       {
 570         nodeOnEdge = findNodeOnEdge( edgeIntPnt, edge );
 571         if ( !nodeOnEdge )
 572         {
 573           if ( shape1.IsNull() && node1 )
 574           {
 575             nodeOnEdge = node1;
 576             node1      = nullptr;
 577             meshDS->MoveNode( nodeOnEdge, edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 578           }
 579           else
 580           {
 581             nodeOnEdge = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 582           }
 583           addNodeOnEdge( nodeOnEdge, edge, paramOnE );
 584         }
 585       }
 586       else
 587       {
 588         _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 589         nodeOnEdge = SMESH_Algo::VertexNode( vertex, meshDS );
 590       }
 591
 592       // create a sub-segment
 593       SMDS_MeshElement* newSeg = _helper.AddEdge( node0, nodeOnEdge );
 594       meshDS->SetMeshElementOnShape( newSeg, shape0 );
 595
 596       SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
 597
 598       node0 = nodeOnEdge;
 599       xyz0  = SMESH_NodeXYZ( node0 );
 600       if ( vertex.IsNull() )
 601         shape0 = edge;
 602       else
 603         shape0 = vertex;
 604     }
 605
 606     else // shape0 is EDGE or VERTEX; look for the next FACE
 607     {    // ------------------------------------------------
 608
 609       if ( !shape1.IsNull() &&
 610            shape1.ShapeType() == TopAbs_FACE &&
 611            _helper.IsSubShape( shape0, shape1 )) // shape0 belongs to FACE shape1
 612       {
 613         SMDS_MeshElement* newSeg = _helper.AddEdge( node0, node1 );
 614         SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
 615         meshDS->SetMeshElementOnShape( newSeg, shape1 );
 616         break;
 617       }
 618       // FACE search
 619       TopoDS_Face face;
 620       double shift = 10 * _tol;
 621       for ( int nbAttemp = 0; face.IsNull() && nbAttemp < 10; ++nbAttemp )
 622       {
 623         xyz0  += segDir * shift;
 624         shift *= 2;
 625         BLSURFPlugin_BLSURF::getProjectionPoint( face, xyz0 );
 626       }
 627       if ( !face.IsNull() )
 628       {
 629         if ( _helper.IsSubShape( shape0, face ))
 630           _faces.Add( face );
 631         else
 632           break;
 633         shape0 = face;
 634       }
 635       else
 636       {
 637         break;
 638       }
 639     }
 640     continue;
 641   } //while ( !shape0.IsSame( shape1 ))
 642
 643   meshDS->RemoveFreeElement( segment, /*submesh=*/nullptr );
 644 }
 645
 646 //================================================================================
 647 /*!
 648  * \brief Find intersection of FACE EDGEs and a segment
 649  *  \param [in] theFace - the FACE
 650  *  \param [in] thePnt0 - first end of the segment
 651  *  \param [in] thePnt1 - last end of the segment
 652  *  \param [out] theFounfEdge - return the intersected EDGE
 653  *  \param [out] theParamOnEdge - return parameter of intersection point on EDGE
 654  *  \return gp_XYZ - point on an EDGE closest to the segment
 655  */
 656 //================================================================================
 657
 658 gp_Pnt BLSURFPlugin_EnforcedMesh1D::getEdgeIntersection( const TopoDS_Shape& theFaceOrEdge,
 659                                                          const gp_XYZ&       thePnt0,
 660                                                          const gp_XYZ&       thePnt1,
 661                                                          TopoDS_Edge &       theFounfEdge,
 662                                                          double &            theParamOnEdge)
 663 {
 664   const double segLen = ( thePnt1 - thePnt0 ).Modulus();
 665   const double maxSegDist2 = segLen * segLen * 0.5 * 0.5;
 666
 667   Handle(Geom_Line) segLine = new Geom_Line( thePnt0, thePnt1 - thePnt0 );
 668   GeomAdaptor_Curve segLineAdpt( segLine, 0, segLen );
 669
 670   TopTools_MapOfShape edges;
 671   double minParamOnSeg = segLen;
 672   gp_Pnt foundPnt;
 673   for ( TopExp_Explorer edgeExp( theFaceOrEdge, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
 674   {
 675     const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
 676     if ( !edges.Add( edge ))
 677       continue;
 678
 679     Extrema_ExtCC extrema( segLineAdpt, BRepAdaptor_Curve( edge ), _tol, _tol );
 680
 681     if ( extrema.IsDone() && !extrema.IsParallel() )
 682       for ( int i = 1, nb = extrema.NbExt(); i <= nb; ++i )
 683         if ( extrema.SquareDistance( i ) < maxSegDist2 )
 684         {
 685           Extrema_POnCurv pOnSeg, pOnEdge;
 686           extrema.Points( i, pOnSeg, pOnEdge );
 687           double paramOnSeg = pOnSeg.Parameter();
 688           if ( 0 < paramOnSeg && paramOnSeg < minParamOnSeg )
 689           {
 690             minParamOnSeg = paramOnSeg;
 691             foundPnt = pOnEdge.Value();
 692             theFounfEdge = edge;
 693             theParamOnEdge = pOnEdge.Parameter();
 694           }
 695         }
 696   }
 697   return foundPnt;
 698 }
 699
 700 //================================================================================
 701 /*!
 702  * \brief Split self-intersecting segments on a given FACE
 703  */
 704 //================================================================================
 705
 706 void BLSURFPlugin_EnforcedMesh1D::splitSelfIntersectingSegments( const TopoDS_Shape & theFace )
 707 {
 708   const TopoDS_Face& face = TopoDS::Face( theFace );
 709
 710   SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
 711   SMESHDS_SubMesh*  sm = meshDS->MeshElements( face );
 712   if ( !sm || sm->NbElements() <= 1 )
 713     return;
 714
 715   // get ordered nodes and segments on the face
 716   SMESH_MeshAlgos::TElemGroupVector edgeBranches;
 717   SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
 718   SMESH_MeshAlgos::Get1DBranches( sm->GetElements(), edgeBranches, nodeBranches );
 719
 720   // create element searcher
 721   SMESH_ElementSearcher*                       elemSearcher;
 722   SMESHUtils::Deleter< SMESH_ElementSearcher > elSearchdeleter;
 723   std::vector< const SMDS_MeshElement* >       foundElems;
 724   {
 725     std::vector< SMDS_ElemIteratorPtr > elemItVec;
 726     for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
 727     {
 728       SMDS_ElemIteratorPtr segIt =
 729         boost::make_shared< SMDS_ElementVectorIterator >( braSegs.begin(), braSegs.end() );
 730       elemItVec.push_back( segIt );
 731     }
 732     typedef SMDS_IteratorOnIterators< const SMDS_MeshElement*,
 733                                       std::vector< SMDS_ElemIteratorPtr > > TVecIterator;
 734     SMDS_ElemIteratorPtr segIt = boost::make_shared< TVecIterator >( elemItVec );
 735
 736     elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, segIt, _tol );
 737     elSearchdeleter._obj = elemSearcher;
 738
 739     // force usage of iterators before they die
 740     elemSearcher->FindElementsByPoint( gp_Pnt( 0,0,1e+20), SMDSAbs_Edge, foundElems );
 741   }
 742
 743
 744   // Find intersecting segments
 745
 746   std::map< const SMDS_MeshElement* , std::vector< const SMDS_MeshNode* > > segInternalNodes;
 747   SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
 748
 749   for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
 750   {
 751     braSegs.push_back( braSegs.back() );
 752     const SMDS_MeshElement* prevSeg = nullptr;
 753
 754     for ( size_t i = 0, nb = braSegs.size() - 1;  i < nb;  ++i )
 755     {
 756       const SMDS_MeshElement*     seg = braSegs[ i ];
 757       const SMDS_MeshElement* nextSeg = braSegs[ i + 1 ];
 758       const SMDS_MeshNode*      node0 = seg->GetNode(0);
 759       const SMDS_MeshNode*      node1 = seg->GetNode(1);
 760
 761       gp_XYZ      xyz0 = SMESH_NodeXYZ( node0 );
 762       gp_XYZ      xyz1 = SMESH_NodeXYZ( node1 );
 763       gp_XYZ middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
 764       double    segLen = ( xyz0 - xyz1 ).Modulus();
 765
 766       foundElems.clear();
 767       elemSearcher->GetElementsInSphere( middlePnt, 0.5 * segLen + _tol, SMDSAbs_Edge, foundElems );
 768
 769       for ( const SMDS_MeshElement* closeSeg : foundElems )
 770       {
 771         if ( closeSeg == prevSeg ||
 772              closeSeg >= seg     ||
 773              closeSeg == nextSeg )
 774           continue;
 775
 776         gp_Pnt   intPnt;
 777         gp_Pnt2d uv;
 778         if ( intersectSegments( seg, closeSeg, face, intPnt, uv ))
 779         {
 780           int i0, i1;
 781           const SMDS_MeshNode* intNode0 = findNode( intPnt, _tol, seg,      segInternalNodes, i0 );
 782           const SMDS_MeshNode* intNode1 = findNode( intPnt, _tol, closeSeg, segInternalNodes, i1 );
 783           if      ( !intNode0 && intNode1 )
 784           {
 785             segInternalNodes[ seg ].push_back( intNode1 );
 786           }
 787           else if ( !intNode1 && intNode0 )
 788           {
 789             segInternalNodes[ closeSeg ].push_back( intNode0 );
 790           }
 791           else if ( intNode1 && intNode0 )
 792           {
 793             if ( intNode1 == intNode0 )
 794               continue;
 795             if ( i0 < 0 && i1 < 0 )
 796             {
 797               nodeGroupsToMerge.push_back  // merge end nodes
 798                 ( std::list< const SMDS_MeshNode* >({ intNode0, intNode1 }));
 799             }
 800             else if ( i0 < 0 )
 801             {
 802               segInternalNodes[ closeSeg ][ i1 ] = intNode0;
 803             }
 804             else if ( i1 < 0 )
 805             {
 806               segInternalNodes[ seg ][ i0 ] = intNode1;
 807             }
 808             else // two internal nodes coincide
 809             {
 810               segInternalNodes[ seg ][ i0 ] = intNode1;
 811               nodeGroupsToMerge.push_back
 812                 ( std::list< const SMDS_MeshNode* >({ intNode1, intNode0 }));
 813             }
 814           }
 815           else // ( !intNode1 && !intNode0 )
 816           {
 817             intNode0 = _helper.AddNode( intPnt.X(), intPnt.Y(), intPnt.Z() );
 818             meshDS->SetNodeOnFace( intNode0, face, uv.X(), uv.Y() );
 819             segInternalNodes[ seg      ].push_back( intNode0 );
 820             segInternalNodes[ closeSeg ].push_back( intNode0 );
 821           }
 822         }
 823       }
 824
 825       prevSeg = seg;
 826
 827     } // loop on segments of a branch
 828   } // loop on branches
 829
 830
 831   findIntersectionWithSeamEdge( face, segInternalNodes ); // on periodic FACE
 832
 833
 834   // Split segments
 835
 836   for ( auto& seg2nodes : segInternalNodes )
 837   {
 838     const SMDS_MeshElement*                 seg = seg2nodes.first;
 839     std::vector< const SMDS_MeshNode* > & nodes = seg2nodes.second;
 840     if ( nodes.empty() ) continue;
 841
 842     const SMDS_MeshNode* n0 = seg->GetNode( 0 );
 843     const SMDS_MeshNode* n1 = seg->GetNode( 1 );
 844     nodes.push_back( n1 );
 845
 846     // sort nodes on the segment
 847     gp_Pnt p0 = SMESH_NodeXYZ( n0 );
 848     std::map< double, const SMDS_MeshNode* > sortedNodes;
 849     for ( SMESH_NodeXYZ pn : nodes )
 850       sortedNodes.insert({ p0.SquareDistance( pn ), pn.Node() });
 851
 852     // make new segments
 853     for ( auto & d2n : sortedNodes )
 854     {
 855       n1 = d2n.second;
 856       SMDS_MeshElement* newSeg = _helper.AddEdge( n0, n1 );
 857       n0 = n1;
 858       meshDS->SetMeshElementOnShape( newSeg, face );
 859       SMESH_MeshEditor::AddToSameGroups( newSeg, seg, meshDS );
 860     }
 861     meshDS->RemoveFreeElement( seg, /*submesh=*/nullptr );
 862   }
 863
 864
 865   // merge equal nodes
 866   SMESH_MeshEditor( _mesh ).MergeNodes( nodeGroupsToMerge );
 867 }
 868
 869 //================================================================================
 870 /*!
 871  * \brief Find intersections of segments with a seam EDGE on a periodic FACE
 872  */
 873 //================================================================================
 874
 875 void BLSURFPlugin_EnforcedMesh1D::findIntersectionWithSeamEdge( const TopoDS_Face & face,
 876                                                                 TNodesOnSeg & segInternalNodes )
 877 {
 878   _helper.SetSubShape( face );
 879   if ( !_helper.HasSeam() )
 880     return;
 881
 882   SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
 883   SMESHDS_SubMesh*  sm = meshDS->MeshElements( face );
 884   if ( !sm || sm->NbElements() == 0 )
 885     return;
 886
 887   TopTools_MapOfShape treatedEdges;
 888   for ( TopExp_Explorer edgeExp( face, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
 889   {
 890     const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
 891     if ( !_helper.IsSeamShape( edge ) || !treatedEdges.Add( edge ))
 892       continue;
 893
 894     for ( SMDS_ElemIteratorPtr setIt = sm->GetElements(); setIt->more(); )
 895     {
 896       const SMDS_MeshElement* segment = setIt->next();
 897       const SMESH_NodeXYZ         pn0 = segment->GetNode( 0 );
 898       const SMESH_NodeXYZ         pn1 = segment->GetNode( 1 );
 899
 900       gp_XY uv0 = _helper.GetNodeUV( face, pn0.Node() );
 901       gp_XY uv1 = _helper.GetNodeUV( face, pn1.Node() );
 902
 903       for ( int iCoo = 1; iCoo <= 2; ++iCoo )
 904         if ( iCoo & _helper.GetPeriodicIndex() )
 905         {
 906           double distParam = Abs( uv0.Coord( iCoo ) - uv1.Coord( iCoo ));
 907           if ( distParam > 0.5 * _helper.GetPeriod( iCoo ))
 908           {
 909             double paramOnE; TopoDS_Edge edge2;
 910             gp_Pnt edgeIntPnt = getEdgeIntersection( edge, pn0, pn1, edge2, paramOnE );
 911             if ( edge2.IsNull() )
 912               continue;
 913
 914             const SMDS_MeshNode* nodeOnSeam = findNodeOnEdge( edgeIntPnt, edge );
 915             bool isOnEdge = nodeOnSeam;
 916
 917             int indexOnSegment = 3;
 918             if ( !nodeOnSeam )
 919               nodeOnSeam = findNode( edgeIntPnt, _tol, segment, segInternalNodes, indexOnSegment );
 920
 921             if ( !nodeOnSeam )
 922             {
 923               nodeOnSeam = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 924             }
 925
 926             if ( !isOnEdge )
 927             {
 928               addNodeOnEdge( nodeOnSeam, edge, paramOnE );
 929             }
 930             if ( indexOnSegment == 3 )
 931               segInternalNodes[ segment ].push_back( nodeOnSeam );
 932           }
 933         }
 934     }
 935   }
 936 }
 937
 938 //================================================================================
 939 /*!
 940  * \brief Intersect two segments
 941  *  \param [in] seg1 - segment 1
 942  *  \param [in] seg2 - segment 2
 943  *  \param [in] face - the FACE on which segments lie
 944  *  \param [out] intPnt - intersection point
 945  *  \param [out] intUV - UV of intersection point on the FACE
 946  *  \return bool - true if intersection found
 947  */
 948 //================================================================================
 949
 950 bool BLSURFPlugin_EnforcedMesh1D::intersectSegments( const SMDS_MeshElement* seg1,
 951                                                      const SMDS_MeshElement* seg2,
 952                                                      const TopoDS_Face&      face,
 953                                                      gp_Pnt &                intPnt,
 954                                                      gp_Pnt2d &              intUV ) const
 955 {
 956   SMESH_NodeXYZ n10 = seg1->GetNode(0);
 957   SMESH_NodeXYZ n11 = seg1->GetNode(1);
 958   SMESH_NodeXYZ n20 = seg2->GetNode(0);
 959   SMESH_NodeXYZ n21 = seg2->GetNode(1);
 960   if ( n10 == n20 || n10 == n21 || n11 == n20 || n10 == n21 )
 961     return false;
 962
 963   gp_Lin lin1( n10, n11 - n10 );
 964   gp_Lin lin2( n20, n21 - n20 );
 965
 966   Extrema_ExtElC extrema( lin1, lin2, Precision::Angular() );
 967
 968   if ( !extrema.IsDone() || extrema.IsParallel() )
 969     return false;
 970
 971   Extrema_POnCurv poc1, poc2;
 972   extrema.Points( 1, poc1, poc2 );
 973
 974   double len1 = lin1.Direction().XYZ() * ( n11 - n10 );
 975   double len2 = lin2.Direction().XYZ() * ( n21 - n20 );
 976   double   u1 = poc1.Parameter();
 977   double   u2 = poc2.Parameter();
 978
 979   if ( u1 < -_tol || u1 > len1 + _tol )
 980     return false;
 981   if ( u2 < -_tol || u2 > len2 + _tol )
 982     return false;
 983
 984   intPnt = 0.5 * ( poc1.Value().XYZ() + poc2.Value().XYZ() );
 985
 986   // compute approximate UV
 987
 988   u1 /= len1;
 989   u2 /= len2;
 990
 991   gp_XY uv10 = _helper.GetNodeUV( face, n10.Node(), n11.Node() );
 992   gp_XY uv11 = _helper.GetNodeUV( face, n11.Node(), n10.Node() );
 993   gp_XY uv1  = uv10 * ( 1 - u1 ) + uv11 * u1;
 994
 995   gp_XY uv20 = _helper.GetNodeUV( face, n20.Node(), n21.Node() );
 996   gp_XY uv21 = _helper.GetNodeUV( face, n21.Node(), n20.Node() );
 997   gp_XY uv2  = uv20 * ( 1 - u2 ) + uv21 * u2;
 998
 999   intUV = 0.5 * ( uv1 + uv2 );
1000
1001   // compute precise UV and XYZ by projecting intPnt to the FACE
1002
1003   Handle(ShapeAnalysis_Surface) surface = _helper.GetSurface( face );
1004   intUV = surface->NextValueOfUV( intUV, intPnt, _tol );
1005   if ( surface->Gap() > Min( len1, len2 ))
1006     intUV = surface->ValueOfUV( intPnt, _tol );
1007
1008   intPnt = surface->Value( intUV );
1009
1010   return true;
1011 }
1012
1013 //================================================================================
1014 /*!
1015  * \brief Setup predicates to detect nodes on FACE boundary
1016  *  \param [in] shape - shape containing FACEs to mesh
1017  */
1018 //================================================================================
1019
1020 void BLSURFPlugin_EnforcedMesh1D::setupPredicates( const TopoDS_Shape& theShape )
1021 {
1022   if ( _onVertexPredicate )
1023     return;
1024
1025   _onVertexPredicate.reset( new TPredicate() );
1026   _onVertexPredicate->SetTolerance( _tol );
1027   _onVertexPredicate->SetMesh( _helper.GetMeshDS() );
1028   {
1029     TopTools_IndexedMapOfShape vertices;
1030     TopExp::MapShapes( theShape, TopAbs_VERTEX, vertices );
1031     TopoDS_Compound vCompound;
1032     BRep_Builder    builder;
1033     builder.MakeCompound( vCompound );
1034     for ( const TopoDS_Shape& v : vertices )
1035       builder.Add( vCompound, v );
1036
1037     _onVertexPredicate->SetShape( vCompound, SMDSAbs_Node );
1038   }
1039
1040
1041   _onEdgePredicate.reset( new TPredicate() );
1042   _onEdgePredicate->SetTolerance( _tol );
1043   _onEdgePredicate->SetMesh( _helper.GetMeshDS() );
1044   {
1045     TopTools_IndexedMapOfShape edges;
1046     TopExp::MapShapes( theShape, TopAbs_EDGE, edges );
1047     TopoDS_Compound eCompound;
1048     BRep_Builder    builder;
1049     builder.MakeCompound( eCompound );
1050     for ( const TopoDS_Shape& e : edges )
1051       builder.Add( eCompound, e );
1052
1053     _onEdgePredicate->SetShape( eCompound, SMDSAbs_Node );
1054   }
1055   return;
1056 }
1057
1058 //================================================================================
1059 /*!
1060  * \brief Find or create a group of edges with given name
1061  */
1062 //================================================================================
1063
1064 SMDS_MeshGroup* BLSURFPlugin_EnforcedMesh1D::getGroup( const std::string& groupName )
1065 {
1066   SMDS_MeshGroup* group = nullptr;
1067   if ( !groupName.empty() )
1068   {
1069     if ( _name2Group.count( groupName ))
1070       return _name2Group[ groupName ];
1071
1072     // find existing group
1073     for ( SMESH_Mesh::GroupIteratorPtr grIt = _mesh->GetGroups(); grIt->more(); )
1074     {
1075       SMESH_Group*     grp = grIt->next();
1076       SMESHDS_Group* grpDS = dynamic_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1077       if ( grpDS &&
1078            grpDS->GetType() == SMDSAbs_Edge &&
1079            groupName == grp->GetName() )
1080       {
1081         _name2Group[ groupName ] = & grpDS->SMDSGroup();
1082         return & grpDS->SMDSGroup();
1083       }
1084     }
1085
1086     // create a new group
1087     SMESH_Group*     grp = _mesh->AddGroup( SMDSAbs_Edge, groupName.c_str() );
1088     SMESHDS_Group* grpDS = static_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1089
1090     group = & grpDS->SMDSGroup();
1091     _name2Group[ groupName ] = group;
1092   }
1093   return group;
1094 }
1095
1096 //================================================================================
1097 /*!
1098  * \brief Look for a node dividing a given EDGE
1099  */
1100 //================================================================================
1101
1102 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::findNodeOnEdge( const gp_Pnt&      p,
1103                                                                   const TopoDS_Edge& edge )
1104 {
1105   // look for an equal node on the EDGE
1106   if ( std::vector< const SMDS_MeshNode* >* nodesOnE = _nodesOnEdge.ChangeSeek( edge ))
1107     for ( const SMDS_MeshNode* n : *nodesOnE )
1108       if ( p.SquareDistance( SMESH_NodeXYZ( n )) < _tol * _tol )
1109       {
1110         return n;
1111       }
1112
1113   return nullptr;
1114 }
1115
1116 //================================================================================
1117 /*!
1118  * \brief Add a node to an EDGE
1119  */
1120 //================================================================================
1121
1122 void BLSURFPlugin_EnforcedMesh1D::addNodeOnEdge( const SMDS_MeshNode* node,
1123                                                  const TopoDS_Edge&   edge,
1124                                                  const double         u)
1125 {
1126   _mesh->GetMeshDS()->SetNodeOnEdge( node, edge, u );
1127
1128   std::vector< const SMDS_MeshNode* > * nodesOnE = _nodesOnEdge.ChangeSeek( edge );
1129   if ( !nodesOnE )
1130     nodesOnE = _nodesOnEdge.Bound( edge, std::vector< const SMDS_MeshNode* >() );
1131
1132   nodesOnE->push_back( node );
1133 }
1134
1135 //================================================================================
1136 /*!
1137  * \brief Create EDGEs by dividing a given EDGE by nodes on it
1138  *  \param [in] edge - the EDGE to divide
1139  *  \param [in] nodes - the nodes to divide by
1140  */
1141 //================================================================================
1142
1143 void BLSURFPlugin_EnforcedMesh1D::
1144
1145 splitEdgeByNodes( TopoDS_Edge edge, const std::vector< const SMDS_MeshNode* >& nodes )
1146 {
1147   if ( nodes.empty() )
1148     return;
1149
1150   edge.Orientation( TopAbs_FORWARD );
1151
1152   TopoDS_Vertex v0 = _helper.IthVertex( 0, edge );
1153   TopoDS_Vertex v1 = _helper.IthVertex( 1, edge );
1154
1155   // create VERTEXes and sort them along the EDGE
1156
1157   std::map< double, TopoDS_Vertex > sortedVertices;
1158
1159   BRepAdaptor_Curve curve( edge );
1160   for ( SMESH_NodeXYZ pn : nodes )
1161   {
1162     gp_Pnt projPnt;
1163     double u;
1164     ShapeAnalysis_Curve().Project( curve, pn, _tol, projPnt, u, false );
1165     projPnt = curve.Value( u );
1166
1167     TopoDS_Vertex v = BRepBuilderAPI_MakeVertex( projPnt );
1168
1169     sortedVertices.insert({ u, v });
1170
1171     _nOnE2Vertex[ pn.Node() ] = v;
1172   }
1173   sortedVertices.insert({ BRep_Tool::Parameter( v1, edge ), v1 });
1174
1175
1176   // create EDGEs
1177
1178   BRep_Builder builder;
1179   std::vector< TopoDS_Edge >& newEdges = *_edgeSplitsOfEdge.Bound( edge, TEdge2Edges::value_type());
1180
1181   double u0 = BRep_Tool::Parameter( v0, edge );
1182   for ( auto& u2v : sortedVertices )
1183   {
1184     double u1 = u2v.first;
1185     v1        = u2v.second;
1186
1187     TopoDS_Shape newShape = edge.EmptyCopied();
1188     TopoDS_Edge  newEdge  = TopoDS::Edge( newShape );
1189     builder.Add( newEdge, v0 );
1190     builder.Add( newEdge, v1 );
1191     builder.Range( newEdge, u0, u1 );
1192     newEdges.push_back( newEdge );
1193
1194     v0 = v1;
1195     u0 = u1;
1196   }
1197
1198   return;
1199 }
1200
1201 //================================================================================
1202 /*!
1203  * \brief  Return true if there are enforced segments on a FACE. Start iteration on segments
1204  */
1205 //================================================================================
1206
1207 bool BLSURFPlugin_EnforcedMesh1D::HasSegmentsOnFace( const TopoDS_Face& face )
1208 {
1209   _segIterator.reset();
1210
1211   if ( _faces.Contains( face ))
1212   {
1213     _currentFace = face;
1214     _segIterator = _helper.GetMeshDS()->MeshElements( face )->GetElements();
1215
1216     _helper.SetSubShape( face );
1217
1218     return _segIterator->more();
1219   }
1220   return false;
1221 }
1222
1223
1224 //================================================================================
1225 /*!
1226  * \brief Return next segment on the FACE
1227  */
1228 //================================================================================
1229
1230 bool BLSURFPlugin_EnforcedMesh1D::NextSegment( Segmemnt &                   seg,
1231                                                TopTools_IndexedMapOfShape & vertexTags )
1232 {
1233   if ( _segIterator && _segIterator->more() )
1234   {
1235     const SMDS_MeshElement* segment = _segIterator->next();
1236
1237     while ( segment->GetType() != SMDSAbs_Edge && _segIterator->more() )
1238       segment = _segIterator->next();
1239     if ( segment->GetType() != SMDSAbs_Edge )
1240       return false;
1241
1242     const SMDS_MeshNode * node[2] = { segment->GetNode(0),
1243                                       segment->GetNode(1) };
1244
1245     seg._tag = _segTag++;
1246
1247     seg._xyz[0] = SMESH_NodeXYZ( node[0]);
1248     seg._xyz[1] = SMESH_NodeXYZ( node[1]);
1249
1250     seg._uv[0] = _helper.GetNodeUV( _currentFace, node[0], node[1] );
1251     seg._uv[1] = _helper.GetNodeUV( _currentFace, node[1], node[0] );
1252
1253     seg._pcurve = new Geom2d_Line( seg._uv[0], ( seg._uv[1] - seg._uv[0] ));
1254
1255     seg._u[0] = 0.0;
1256     seg._u[1] = ( seg._uv[1] - seg._uv[0] ).Modulus();
1257
1258     for ( int i = 0; i < 2; ++i )
1259     {
1260       auto n2v = _nOnE2Vertex.find( node[i] ); // find VERTEX by node
1261
1262       if ( n2v != _nOnE2Vertex.end() )
1263         seg._vTag[i] = vertexTags.Add( n2v->second );
1264       else
1265         seg._vTag[i] = _nodeTag0 + _nodeTags.Add( node[i] );
1266     }
1267
1268     if ( !_isQuadratic )
1269       _mesh->GetMeshDS()->UnSetMeshElementOnShape( segment, _currentFace );
1270
1271     return true;
1272   }
1273   return false;
1274 }
1275
1276 //================================================================================
1277 /*!
1278  * \brief Return enforced node by the tag that was returned by Segmemnt::_vTag[i]
1279  */
1280 //================================================================================
1281
1282 const SMDS_MeshNode*
1283 BLSURFPlugin_EnforcedMesh1D::GetNodeByTag( int                                tag,
1284                                            const TopTools_IndexedMapOfShape & vertexTags )
1285 {
1286   const SMDS_MeshNode* node = nullptr;
1287
1288   if ( tag <= vertexTags.Size() && !_nOnE2Vertex.empty() )
1289   {
1290     const TopoDS_Shape& vertex = vertexTags( tag );
1291
1292     for ( auto n2v = _nOnE2Vertex.begin(); n2v != _nOnE2Vertex.end(); ++n2v )
1293       if ( vertex.IsSame( n2v->second ))
1294       {
1295         node = n2v->first;
1296         _nOnE2Vertex.erase( n2v );
1297         return node;
1298       }
1299   }
1300
1301   tag -= _nodeTag0;
1302
1303   bool isValid = ( 0 < tag && tag <= _nodeTags.Size() );
1304   if ( isValid )
1305     node = _nodeTags( tag );
1306
1307   return node;
1308 }
1309
1310 //================================================================================
1311 /*!
1312  * \brief Return true if a tag corresponds to the tag of enforced segment
1313  *  that was returned by Segment::_tag
1314  */
1315 //================================================================================
1316
1317 bool BLSURFPlugin_EnforcedMesh1D::IsSegmentTag( int tag ) const
1318 {
1319   return ( _segTag0 <= tag && tag <= _segTag );
1320 }
1321
1322 //================================================================================
1323 /*!
1324  * \brief Return tag of EDGE by tags of its splits
1325  */
1326 //================================================================================
1327
1328 int BLSURFPlugin_EnforcedMesh1D::GetTagOfSplitEdge( int splitTag ) const
1329 {
1330   auto sTag2eTag = _splitTags2EdgeTag.find( splitTag );
1331   return ( sTag2eTag == _splitTags2EdgeTag.end() ) ? splitTag : sTag2eTag->second;
1332 }