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 Add a vertex on EDGE
 244  */
 245 //================================================================================
 246
 247 void BLSURFPlugin_EnforcedMesh1D::AddVertexOnEdge( const double* theXYZ )
 248 {
 249   // setup predicates to find the supporting EDGE
 250   setupPredicates( _shape );
 251
 252   SMESHDS_Mesh*         meshDS = _mesh->GetMeshDS();
 253   const SMDS_MeshNode* nodeOnE = meshDS->AddNode( theXYZ[0], theXYZ[1], theXYZ[2] );
 254
 255   // check if enfNode is on VERTEX
 256   bool toRemove = true;
 257   TopoDS_Vertex vertex;
 258   TopoDS_Edge   edge;
 259   if ( _onVertexPredicate->IsSatisfy( nodeOnE, &vertex ))
 260   {
 261     toRemove = SMESH_Algo::VertexNode( vertex, meshDS );
 262     if ( !toRemove )
 263       meshDS->SetNodeOnVertex( nodeOnE, vertex );
 264   }
 265   // find the EDGE supporting theXYZ
 266   else if ( _onEdgePredicate->IsSatisfy( nodeOnE, &edge ))
 267   {
 268     gp_Pnt pnt( theXYZ[0], theXYZ[1], theXYZ[2] );
 269     toRemove = findNodeOnEdge( pnt, edge );
 270     if ( !toRemove )
 271       addNodeOnEdge( nodeOnE, edge );
 272   }
 273
 274   if ( toRemove )
 275     meshDS->RemoveFreeNode( nodeOnE, /*submesh=*/nullptr, /*fromGroup=*/false );
 276 }
 277
 278 //================================================================================
 279 /*!
 280  * \brief Return EDGEs resulted from division of FACE boundary by enforced segments
 281  *        and enforced vertices
 282  *  \param [in] edge - possibly divided EDGE
 283  *  \param [out] splits - split EDGEs
 284  *  \return bool - true if the EDGE is split
 285  */
 286 //================================================================================
 287
 288 bool BLSURFPlugin_EnforcedMesh1D::GetSplitsOfEdge( const TopoDS_Edge&           edge,
 289                                                    std::vector< TopoDS_Edge > & splits,
 290                                                    TopTools_IndexedMapOfShape & edgeTags )
 291 {
 292   if ( _nodesOnEdge.IsBound( edge )) // divide the EDGE
 293   {
 294     splitEdgeByNodes( edge, _nodesOnEdge( edge ));
 295
 296     _nodesOnEdge.UnBind( edge );
 297   }
 298
 299   // return splits
 300
 301   std::vector< TopoDS_Edge > * splitsInMap = _edgeSplitsOfEdge.ChangeSeek( edge );
 302   if ( !splitsInMap )
 303     return false;
 304
 305   splits.insert( splits.end(), splitsInMap->begin(), splitsInMap->end() );
 306
 307   int eTag = edgeTags.Add( edge );
 308
 309   size_t index = splits.size() - 1;
 310   for ( size_t i = 0; i < splitsInMap->size(); ++i, --index )
 311   {
 312     int splitTag = edgeTags.Add( splits[ index ]);
 313     _splitTags2EdgeTag[ splitTag ] = eTag;
 314
 315     if ( edge.Orientation() == TopAbs_REVERSED )
 316       splits[ index ].Reverse();
 317   }
 318
 319   return true;
 320 }
 321
 322
 323 //================================================================================
 324 /*!
 325  * \brief Add 1D elements to the mesh
 326  */
 327 //================================================================================
 328
 329 void BLSURFPlugin_EnforcedMesh1D::
 330
 331 copyEnforcedMesh( const BLSURFPlugin_Hypothesis::EnforcedMesh& theEnfMesh,
 332                   const BLSURFPlugin_Hypothesis*               theHyp,
 333                   const TopoDS_Shape&                          theShape)
 334 {
 335   SMESH_Mesh* mesh1D;
 336   SMDS_ElemIteratorPtr segIt = theHyp->GetEnforcedSegments( theEnfMesh, mesh1D );
 337   if ( !segIt->more() )
 338     return;
 339
 340   // setup predicates to detect nodes on FACE boundary
 341   setupPredicates( theShape );
 342
 343   SMDS_MeshGroup* group = getGroup( theEnfMesh._groupName );
 344   SMESHDS_Mesh*  meshDS = _helper.GetMeshDS();
 345
 346   // get ordered nodes and segments of theEnfMesh
 347   SMESH_MeshAlgos::TElemGroupVector edgeBranches;
 348   SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
 349   SMESH_MeshAlgos::Get1DBranches( segIt, edgeBranches, nodeBranches );
 350
 351
 352   // Copy nodes and segments from an enforced mesh to my mesh
 353
 354   TopoDS_Shape vertex, edge;
 355
 356   // first treat ends of branches that can be shared by branches
 357   for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
 358   {
 359     std::vector< const SMDS_MeshNode* > &   braNodes = nodeBranches[ iB ];
 360     std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
 361
 362     for ( int isLast = 0; isLast < 2; ++isLast )
 363     {
 364       const SMDS_MeshNode* enfNode = isLast ? braNodes.back() : braNodes[0];
 365       if ( meshDS->Contains( enfNode ))
 366         continue; // already in my mesh
 367
 368       const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
 369       if ( !newNode )
 370         orientSegments( braSegs, braNodes, isLast ? 0 : braNodes.size() - 1, mesh1D );
 371
 372       // replace enfNode at branch ends by newNode
 373       SMESH_NodeXYZ enfPnt( newNode ? newNode : enfNode );
 374       for ( std::vector< const SMDS_MeshNode* > & braNodes : nodeBranches )
 375       {
 376         for ( int isLast = 0; isLast < 2; ++isLast )
 377         {
 378           const SMDS_MeshNode* & endNode = isLast ? braNodes.back() : braNodes[0];
 379           if ( endNode == enfNode || enfPnt.SquareDistance( endNode ) < _tol*_tol )
 380             endNode = newNode;
 381         }
 382       }
 383       continue;
 384     }  // loop on branch ends
 385   } // loop on nodeBranches
 386
 387   // copy nodes and segments
 388
 389   for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
 390   {
 391     std::vector< const SMDS_MeshNode* > &   braNodes = nodeBranches[ iB ];
 392     std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
 393
 394     // copy nodes of the branch
 395     for ( size_t i = 0; i < braNodes.size(); ++i )
 396     {
 397       const SMDS_MeshNode* & enfNode = braNodes[ i ];
 398       const SMDS_MeshNode*   newNode = copyEnforcedNode( enfNode );
 399
 400       if ( !newNode ) // orient segments to be able to get enforced not projected node
 401         orientSegments( braSegs, braNodes, i, mesh1D );
 402
 403       enfNode = newNode;
 404     }
 405
 406     // copy segments of the branch
 407     for ( size_t i = 0; i < braSegs.size(); ++i )
 408     {
 409       //braSegs[ i ] = nullptr;
 410
 411       const SMDS_MeshNode* node0 = braNodes[ i     ];
 412       const SMDS_MeshNode* node1 = braNodes[ i + 1 ];
 413       if ( !node0 && !node1 )
 414         continue;
 415
 416       TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
 417       TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
 418       if ( shape0.IsNull() && shape1.IsNull() )
 419         continue;
 420
 421       if ( !node0 && shape1.ShapeType() != TopAbs_FACE )
 422         continue;
 423       if ( !node1 && shape0.ShapeType() != TopAbs_FACE )
 424         continue;
 425
 426       if ( !node0 || !node1 ) // create missing node at location of enforced node projected nowhere
 427       {
 428         SMESH_NodeXYZ  pn = braSegs[i]->GetNode( !node1 );
 429         ( node0 ? node1 : node0 ) = _helper.AddNode( pn->X(), pn->Y(), pn->Z() );
 430       }
 431
 432       SMDS_MeshEdge* newSeg = _helper.AddEdge( node0, node1 );
 433       if ( group )
 434         group->Add( newSeg );
 435       braSegs[ i ] = newSeg;
 436
 437       // check if the both nodes are on the same FACE
 438       TopoDS_Shape face = shape0;
 439       if ( !shape0.IsSame( shape1 ) && !shape0.IsNull() && !shape1.IsNull() )
 440       {
 441         if ( shape0.ShapeType() == TopAbs_FACE &&
 442              _helper.IsSubShape( shape1, shape0 ))
 443         {
 444           face = shape0;
 445         }
 446         else if ( shape1.ShapeType() == TopAbs_FACE &&
 447                   _helper.IsSubShape( shape0, shape1 ))
 448         {
 449           face = shape1;
 450         }
 451         else // try to find a FACE by projecting a segment middle point
 452         {
 453           face.Nullify();
 454           gp_Pnt middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
 455           //BLSURFPlugin_BLSURF::projectionPoint projPnt =
 456           BLSURFPlugin_BLSURF::getProjectionPoint( TopoDS::Face( face ), middlePnt );
 457
 458           if ( !face.IsNull() &&
 459                ( !_helper.IsSubShape( shape0, face ) ||
 460                  !_helper.IsSubShape( shape1, face ) ))
 461             face.Nullify();
 462         }
 463       }
 464       if ( !face.IsNull() && face.ShapeType() == TopAbs_FACE )
 465       {
 466         meshDS->SetMeshElementOnShape( newSeg, face );
 467         _faces.Add( face );
 468       }
 469
 470       if ( face.IsNull() || shape0.IsNull() || shape1.IsNull() )
 471       {
 472         _segmentsOnSeveralShapes.push_back( newSeg );
 473       }
 474
 475     } // loop on branch segments
 476     continue;
 477   } // loop on branches
 478
 479   return;
 480 }
 481
 482 //================================================================================
 483 /*!
 484  * \brief Create a copy of a node of enforced mesh in my mesh
 485  *  \param [in] enfNode - enforced node
 486  *  \return const SMDS_MeshNode* - a node in my mesh
 487  */
 488 //================================================================================
 489
 490 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::copyEnforcedNode( const SMDS_MeshNode* enfNode )
 491 {
 492   SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
 493
 494   if ( !enfNode || meshDS->Contains( enfNode ))
 495     return enfNode; // already in my mesh
 496
 497   SMESH_NodeXYZ enfPnt = enfNode;
 498
 499   const SMDS_MeshNode* newNode = nullptr;
 500
 501   // check if enfNode is on VERTEX
 502   TopoDS_Vertex vertex;
 503   if ( _onVertexPredicate->IsSatisfy( enfNode, &vertex ))
 504   {
 505     _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 506     newNode = SMESH_Algo::VertexNode( vertex, meshDS );
 507   }
 508
 509   // check if enfNode is on EDGE
 510   bool setNodeOnEdge = false;
 511   TopoDS_Edge edge;
 512   if ( !newNode )
 513   {
 514     setNodeOnEdge = _onEdgePredicate->IsSatisfy( enfNode, &edge );
 515     if ( setNodeOnEdge )
 516     {
 517       newNode = findNodeOnEdge( enfPnt, edge );
 518       setNodeOnEdge = !newNode;
 519     }
 520   }
 521
 522   // create a new node and set it on FACE
 523   if ( !newNode )
 524   {
 525     TopoDS_Face face;
 526     BLSURFPlugin_BLSURF::projectionPoint projPnt =
 527       BLSURFPlugin_BLSURF::getProjectionPoint( face, enfPnt, /*allowStateON=*/true );
 528     if ( face.IsNull() ) return newNode;
 529
 530     if ( projPnt.state == TopAbs_ON )
 531     {
 532       SMDS_MeshNode* projNode = const_cast< SMDS_MeshNode* >( enfNode );
 533       projNode->setXYZ( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
 534       vertex.Nullify();
 535       edge.Nullify();
 536       if ( _onVertexPredicate->IsSatisfy( projNode, &vertex ))
 537       {
 538         _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 539         newNode = SMESH_Algo::VertexNode( vertex, meshDS );
 540       }
 541       else if (( setNodeOnEdge = _onEdgePredicate->IsSatisfy( projNode, &edge )))
 542       {
 543         newNode = findNodeOnEdge( projPnt.xyz, edge );
 544         setNodeOnEdge = !newNode;
 545       }
 546       projNode->setXYZ( enfPnt.X(), enfPnt.Y(), enfPnt.Z() );
 547     }
 548
 549     if ( !newNode )
 550       newNode = meshDS->AddNode( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
 551
 552     if ( vertex.IsNull() && edge.IsNull() )
 553       meshDS->SetNodeOnFace( newNode, face, projPnt.uv.X(), projPnt.uv.Y() );
 554   }
 555
 556   // set the new node on EDGE
 557   if ( newNode && setNodeOnEdge )
 558   {
 559     addNodeOnEdge( newNode, edge );
 560   }
 561
 562   return newNode;
 563 }
 564
 565 //================================================================================
 566 /*!
 567  * \brief Split a segment whose nodes are on different FACEs into smaller parts
 568  *        lying each on one FACE
 569  */
 570 //================================================================================
 571
 572 void BLSURFPlugin_EnforcedMesh1D::splitSegmentOnSeveralShapes( const SMDS_MeshElement* segment )
 573 {
 574   const SMDS_MeshNode* node0 = segment->GetNode(0);
 575   const SMDS_MeshNode* node1 = segment->GetNode(1);
 576   SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
 577
 578   TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
 579   TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
 580
 581   if ( shape0.IsNull() )
 582   {
 583     std::swap( node0, node1 );
 584     shape0 = shape1;
 585     shape1.Nullify();
 586   }
 587
 588   gp_XYZ   xyz0 = SMESH_NodeXYZ( node0 );
 589   gp_XYZ   xyz1 = SMESH_NodeXYZ( node1 );
 590   gp_XYZ segDir = ( xyz1 - xyz0 ).Normalized();
 591
 592   //std::map< double, const SMDS_MeshNode* > mediumNodes; // nodes splitting the segment
 593
 594   while ( !shape0.IsSame( shape1 )) // move along the segment till shape1
 595   {
 596     if ( shape0.ShapeType() == TopAbs_FACE ) // make a node on an EDGE of the FACE
 597     {                                        // ----------------------------------
 598       TopoDS_Edge edge;
 599       double paramOnE;
 600       gp_Pnt edgeIntPnt = getEdgeIntersection( shape0, xyz0, xyz1, edge, paramOnE );
 601       if ( edge.IsNull() )
 602         break;
 603
 604       // check if edgeIntPnt in on VERTEX
 605       TopoDS_Vertex vertex;
 606       for ( int iV = 0; iV < 2 &&  vertex.IsNull(); ++iV )
 607       {
 608         TopoDS_Vertex v = _helper.IthVertex( iV, edge );
 609         if ( edgeIntPnt.SquareDistance( BRep_Tool::Pnt( v )) < _tol *_tol )
 610           vertex = v;
 611       }
 612
 613       // make a node on the EDGE
 614       const SMDS_MeshNode* nodeOnEdge = nullptr;
 615       if ( vertex.IsNull() )
 616       {
 617         nodeOnEdge = findNodeOnEdge( edgeIntPnt, edge );
 618         if ( !nodeOnEdge )
 619         {
 620           if ( shape1.IsNull() && node1 )
 621           {
 622             nodeOnEdge = node1;
 623             node1      = nullptr;
 624             meshDS->MoveNode( nodeOnEdge, edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 625           }
 626           else
 627           {
 628             nodeOnEdge = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 629           }
 630           addNodeOnEdge( nodeOnEdge, edge, paramOnE );
 631         }
 632       }
 633       else
 634       {
 635         _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 636         nodeOnEdge = SMESH_Algo::VertexNode( vertex, meshDS );
 637       }
 638
 639       // create a sub-segment
 640       SMDS_MeshElement* newSeg = _helper.AddEdge( node0, nodeOnEdge );
 641       meshDS->SetMeshElementOnShape( newSeg, shape0 );
 642
 643       SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
 644
 645       node0 = nodeOnEdge;
 646       xyz0  = SMESH_NodeXYZ( node0 );
 647       if ( vertex.IsNull() )
 648         shape0 = edge;
 649       else
 650         shape0 = vertex;
 651     }
 652
 653     else // shape0 is EDGE or VERTEX; look for the next FACE
 654     {    // ------------------------------------------------
 655
 656       if ( !shape1.IsNull() &&
 657            shape1.ShapeType() == TopAbs_FACE &&
 658            _helper.IsSubShape( shape0, shape1 )) // shape0 belongs to FACE shape1
 659       {
 660         SMDS_MeshElement* newSeg = _helper.AddEdge( node0, node1 );
 661         SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
 662         meshDS->SetMeshElementOnShape( newSeg, shape1 );
 663         break;
 664       }
 665       // FACE search
 666       TopoDS_Face face;
 667       double shift = 10 * _tol;
 668       for ( int nbAttemp = 0; face.IsNull() && nbAttemp < 10; ++nbAttemp )
 669       {
 670         xyz0  += segDir * shift;
 671         shift *= 2;
 672         BLSURFPlugin_BLSURF::getProjectionPoint( face, xyz0 );
 673       }
 674       if ( !face.IsNull() )
 675       {
 676         if ( _helper.IsSubShape( shape0, face ))
 677           _faces.Add( face );
 678         else
 679           break;
 680         shape0 = face;
 681       }
 682       else
 683       {
 684         break;
 685       }
 686     }
 687     continue;
 688   } //while ( !shape0.IsSame( shape1 ))
 689
 690   meshDS->RemoveFreeElement( segment, /*submesh=*/nullptr );
 691 }
 692
 693 //================================================================================
 694 /*!
 695  * \brief Find intersection of FACE EDGEs and a segment
 696  *  \param [in] theFace - the FACE
 697  *  \param [in] thePnt0 - first end of the segment
 698  *  \param [in] thePnt1 - last end of the segment
 699  *  \param [out] theFounfEdge - return the intersected EDGE
 700  *  \param [out] theParamOnEdge - return parameter of intersection point on EDGE
 701  *  \return gp_XYZ - point on an EDGE closest to the segment
 702  */
 703 //================================================================================
 704
 705 gp_Pnt BLSURFPlugin_EnforcedMesh1D::getEdgeIntersection( const TopoDS_Shape& theFaceOrEdge,
 706                                                          const gp_XYZ&       thePnt0,
 707                                                          const gp_XYZ&       thePnt1,
 708                                                          TopoDS_Edge &       theFounfEdge,
 709                                                          double &            theParamOnEdge)
 710 {
 711   const double segLen = ( thePnt1 - thePnt0 ).Modulus();
 712   const double maxSegDist2 = segLen * segLen * 0.5 * 0.5;
 713
 714   Handle(Geom_Line) segLine = new Geom_Line( thePnt0, thePnt1 - thePnt0 );
 715   GeomAdaptor_Curve segLineAdpt( segLine, 0, segLen );
 716
 717   TopTools_MapOfShape edges;
 718   double minParamOnSeg = segLen;
 719   gp_Pnt foundPnt;
 720   for ( TopExp_Explorer edgeExp( theFaceOrEdge, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
 721   {
 722     const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
 723     if ( !edges.Add( edge ))
 724       continue;
 725
 726     Extrema_ExtCC extrema( segLineAdpt, BRepAdaptor_Curve( edge ), _tol, _tol );
 727
 728     if ( extrema.IsDone() && !extrema.IsParallel() )
 729       for ( int i = 1, nb = extrema.NbExt(); i <= nb; ++i )
 730         if ( extrema.SquareDistance( i ) < maxSegDist2 )
 731         {
 732           Extrema_POnCurv pOnSeg, pOnEdge;
 733           extrema.Points( i, pOnSeg, pOnEdge );
 734           double paramOnSeg = pOnSeg.Parameter();
 735           if ( 0 < paramOnSeg && paramOnSeg < minParamOnSeg )
 736           {
 737             minParamOnSeg = paramOnSeg;
 738             foundPnt = pOnEdge.Value();
 739             theFounfEdge = edge;
 740             theParamOnEdge = pOnEdge.Parameter();
 741           }
 742         }
 743   }
 744   return foundPnt;
 745 }
 746
 747 //================================================================================
 748 /*!
 749  * \brief Split self-intersecting segments on a given FACE
 750  */
 751 //================================================================================
 752
 753 void BLSURFPlugin_EnforcedMesh1D::splitSelfIntersectingSegments( const TopoDS_Shape & theFace )
 754 {
 755   const TopoDS_Face& face = TopoDS::Face( theFace );
 756
 757   SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
 758   SMESHDS_SubMesh*  sm = meshDS->MeshElements( face );
 759   if ( !sm || sm->NbElements() <= 1 )
 760     return;
 761
 762   // get ordered nodes and segments on the face
 763   SMESH_MeshAlgos::TElemGroupVector edgeBranches;
 764   SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
 765   SMESH_MeshAlgos::Get1DBranches( sm->GetElements(), edgeBranches, nodeBranches );
 766
 767   // create element searcher
 768   SMESH_ElementSearcher*                       elemSearcher;
 769   SMESHUtils::Deleter< SMESH_ElementSearcher > elSearchdeleter;
 770   std::vector< const SMDS_MeshElement* >       foundElems;
 771   {
 772     std::vector< SMDS_ElemIteratorPtr > elemItVec;
 773     for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
 774     {
 775       SMDS_ElemIteratorPtr segIt =
 776         boost::make_shared< SMDS_ElementVectorIterator >( braSegs.begin(), braSegs.end() );
 777       elemItVec.push_back( segIt );
 778     }
 779     typedef SMDS_IteratorOnIterators< const SMDS_MeshElement*,
 780                                       std::vector< SMDS_ElemIteratorPtr > > TVecIterator;
 781     SMDS_ElemIteratorPtr segIt = boost::make_shared< TVecIterator >( elemItVec );
 782
 783     elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, segIt, _tol );
 784     elSearchdeleter._obj = elemSearcher;
 785
 786     // force usage of iterators before they die
 787     elemSearcher->FindElementsByPoint( gp_Pnt( 0,0,1e+20), SMDSAbs_Edge, foundElems );
 788   }
 789
 790
 791   // Find intersecting segments
 792
 793   std::map< const SMDS_MeshElement* , std::vector< const SMDS_MeshNode* > > segInternalNodes;
 794   SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
 795
 796   for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
 797   {
 798     braSegs.push_back( braSegs.back() );
 799     const SMDS_MeshElement* prevSeg = nullptr;
 800
 801     for ( size_t i = 0, nb = braSegs.size() - 1;  i < nb;  ++i )
 802     {
 803       const SMDS_MeshElement*     seg = braSegs[ i ];
 804       const SMDS_MeshElement* nextSeg = braSegs[ i + 1 ];
 805       const SMDS_MeshNode*      node0 = seg->GetNode(0);
 806       const SMDS_MeshNode*      node1 = seg->GetNode(1);
 807
 808       gp_XYZ      xyz0 = SMESH_NodeXYZ( node0 );
 809       gp_XYZ      xyz1 = SMESH_NodeXYZ( node1 );
 810       gp_XYZ middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
 811       double    segLen = ( xyz0 - xyz1 ).Modulus();
 812
 813       foundElems.clear();
 814       elemSearcher->GetElementsInSphere( middlePnt, 0.5 * segLen + _tol, SMDSAbs_Edge, foundElems );
 815
 816       for ( const SMDS_MeshElement* closeSeg : foundElems )
 817       {
 818         if ( closeSeg == prevSeg ||
 819              closeSeg >= seg     ||
 820              closeSeg == nextSeg )
 821           continue;
 822
 823         gp_Pnt   intPnt;
 824         gp_Pnt2d uv;
 825         if ( intersectSegments( seg, closeSeg, face, intPnt, uv ))
 826         {
 827           int i0, i1;
 828           const SMDS_MeshNode* intNode0 = findNode( intPnt, _tol, seg,      segInternalNodes, i0 );
 829           const SMDS_MeshNode* intNode1 = findNode( intPnt, _tol, closeSeg, segInternalNodes, i1 );
 830           if      ( !intNode0 && intNode1 )
 831           {
 832             segInternalNodes[ seg ].push_back( intNode1 );
 833           }
 834           else if ( !intNode1 && intNode0 )
 835           {
 836             segInternalNodes[ closeSeg ].push_back( intNode0 );
 837           }
 838           else if ( intNode1 && intNode0 )
 839           {
 840             if ( intNode1 == intNode0 )
 841               continue;
 842             if ( i0 < 0 && i1 < 0 )
 843             {
 844               nodeGroupsToMerge.push_back  // merge end nodes
 845                 ( std::list< const SMDS_MeshNode* >({ intNode0, intNode1 }));
 846             }
 847             else if ( i0 < 0 )
 848             {
 849               segInternalNodes[ closeSeg ][ i1 ] = intNode0;
 850             }
 851             else if ( i1 < 0 )
 852             {
 853               segInternalNodes[ seg ][ i0 ] = intNode1;
 854             }
 855             else // two internal nodes coincide
 856             {
 857               segInternalNodes[ seg ][ i0 ] = intNode1;
 858               nodeGroupsToMerge.push_back
 859                 ( std::list< const SMDS_MeshNode* >({ intNode1, intNode0 }));
 860             }
 861           }
 862           else // ( !intNode1 && !intNode0 )
 863           {
 864             intNode0 = _helper.AddNode( intPnt.X(), intPnt.Y(), intPnt.Z() );
 865             meshDS->SetNodeOnFace( intNode0, face, uv.X(), uv.Y() );
 866             segInternalNodes[ seg      ].push_back( intNode0 );
 867             segInternalNodes[ closeSeg ].push_back( intNode0 );
 868           }
 869         }
 870       }
 871
 872       prevSeg = seg;
 873
 874     } // loop on segments of a branch
 875   } // loop on branches
 876
 877
 878   findIntersectionWithSeamEdge( face, segInternalNodes ); // on periodic FACE
 879
 880
 881   // Split segments
 882
 883   for ( auto& seg2nodes : segInternalNodes )
 884   {
 885     const SMDS_MeshElement*                 seg = seg2nodes.first;
 886     std::vector< const SMDS_MeshNode* > & nodes = seg2nodes.second;
 887     if ( nodes.empty() ) continue;
 888
 889     const SMDS_MeshNode* n0 = seg->GetNode( 0 );
 890     const SMDS_MeshNode* n1 = seg->GetNode( 1 );
 891     nodes.push_back( n1 );
 892
 893     // sort nodes on the segment
 894     gp_Pnt p0 = SMESH_NodeXYZ( n0 );
 895     std::map< double, const SMDS_MeshNode* > sortedNodes;
 896     for ( SMESH_NodeXYZ pn : nodes )
 897       sortedNodes.insert({ p0.SquareDistance( pn ), pn.Node() });
 898
 899     // make new segments
 900     for ( auto & d2n : sortedNodes )
 901     {
 902       n1 = d2n.second;
 903       SMDS_MeshElement* newSeg = _helper.AddEdge( n0, n1 );
 904       n0 = n1;
 905       meshDS->SetMeshElementOnShape( newSeg, face );
 906       SMESH_MeshEditor::AddToSameGroups( newSeg, seg, meshDS );
 907     }
 908     meshDS->RemoveFreeElement( seg, /*submesh=*/nullptr );
 909   }
 910
 911
 912   // merge equal nodes
 913   SMESH_MeshEditor( _mesh ).MergeNodes( nodeGroupsToMerge );
 914 }
 915
 916 //================================================================================
 917 /*!
 918  * \brief Find intersections of segments with a seam EDGE on a periodic FACE
 919  */
 920 //================================================================================
 921
 922 void BLSURFPlugin_EnforcedMesh1D::findIntersectionWithSeamEdge( const TopoDS_Face & face,
 923                                                                 TNodesOnSeg & segInternalNodes )
 924 {
 925   _helper.SetSubShape( face );
 926   if ( !_helper.HasSeam() )
 927     return;
 928
 929   SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
 930   SMESHDS_SubMesh*  sm = meshDS->MeshElements( face );
 931   if ( !sm || sm->NbElements() == 0 )
 932     return;
 933
 934   TopTools_MapOfShape treatedEdges;
 935   for ( TopExp_Explorer edgeExp( face, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
 936   {
 937     const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
 938     if ( !_helper.IsSeamShape( edge ) || !treatedEdges.Add( edge ))
 939       continue;
 940
 941     for ( SMDS_ElemIteratorPtr setIt = sm->GetElements(); setIt->more(); )
 942     {
 943       const SMDS_MeshElement* segment = setIt->next();
 944       const SMESH_NodeXYZ         pn0 = segment->GetNode( 0 );
 945       const SMESH_NodeXYZ         pn1 = segment->GetNode( 1 );
 946
 947       gp_XY uv0 = _helper.GetNodeUV( face, pn0.Node() );
 948       gp_XY uv1 = _helper.GetNodeUV( face, pn1.Node() );
 949
 950       for ( int iCoo = 1; iCoo <= 2; ++iCoo )
 951         if ( iCoo & _helper.GetPeriodicIndex() )
 952         {
 953           double distParam = Abs( uv0.Coord( iCoo ) - uv1.Coord( iCoo ));
 954           if ( distParam > 0.5 * _helper.GetPeriod( iCoo ))
 955           {
 956             double paramOnE; TopoDS_Edge edge2;
 957             gp_Pnt edgeIntPnt = getEdgeIntersection( edge, pn0, pn1, edge2, paramOnE );
 958             if ( edge2.IsNull() )
 959               continue;
 960
 961             const SMDS_MeshNode* nodeOnSeam = findNodeOnEdge( edgeIntPnt, edge );
 962             bool isOnEdge = nodeOnSeam;
 963
 964             int indexOnSegment = 3;
 965             if ( !nodeOnSeam )
 966               nodeOnSeam = findNode( edgeIntPnt, _tol, segment, segInternalNodes, indexOnSegment );
 967
 968             if ( !nodeOnSeam )
 969             {
 970               nodeOnSeam = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
 971             }
 972
 973             if ( !isOnEdge )
 974             {
 975               addNodeOnEdge( nodeOnSeam, edge, paramOnE );
 976             }
 977             if ( indexOnSegment == 3 )
 978               segInternalNodes[ segment ].push_back( nodeOnSeam );
 979           }
 980         }
 981     }
 982   }
 983 }
 984
 985 //================================================================================
 986 /*!
 987  * \brief Intersect two segments
 988  *  \param [in] seg1 - segment 1
 989  *  \param [in] seg2 - segment 2
 990  *  \param [in] face - the FACE on which segments lie
 991  *  \param [out] intPnt - intersection point
 992  *  \param [out] intUV - UV of intersection point on the FACE
 993  *  \return bool - true if intersection found
 994  */
 995 //================================================================================
 996
 997 bool BLSURFPlugin_EnforcedMesh1D::intersectSegments( const SMDS_MeshElement* seg1,
 998                                                      const SMDS_MeshElement* seg2,
 999                                                      const TopoDS_Face&      face,
1000                                                      gp_Pnt &                intPnt,
1001                                                      gp_Pnt2d &              intUV ) const
1002 {
1003   SMESH_NodeXYZ n10 = seg1->GetNode(0);
1004   SMESH_NodeXYZ n11 = seg1->GetNode(1);
1005   SMESH_NodeXYZ n20 = seg2->GetNode(0);
1006   SMESH_NodeXYZ n21 = seg2->GetNode(1);
1007   if ( n10 == n20 || n10 == n21 || n11 == n20 || n10 == n21 )
1008     return false;
1009
1010   gp_Lin lin1( n10, n11 - n10 );
1011   gp_Lin lin2( n20, n21 - n20 );
1012
1013   Extrema_ExtElC extrema( lin1, lin2, Precision::Angular() );
1014
1015   if ( !extrema.IsDone() || extrema.IsParallel() )
1016     return false;
1017
1018   Extrema_POnCurv poc1, poc2;
1019   extrema.Points( 1, poc1, poc2 );
1020
1021   double len1 = lin1.Direction().XYZ() * ( n11 - n10 );
1022   double len2 = lin2.Direction().XYZ() * ( n21 - n20 );
1023   double   u1 = poc1.Parameter();
1024   double   u2 = poc2.Parameter();
1025
1026   if ( u1 < -_tol || u1 > len1 + _tol )
1027     return false;
1028   if ( u2 < -_tol || u2 > len2 + _tol )
1029     return false;
1030
1031   intPnt = 0.5 * ( poc1.Value().XYZ() + poc2.Value().XYZ() );
1032
1033   // compute approximate UV
1034
1035   u1 /= len1;
1036   u2 /= len2;
1037
1038   gp_XY uv10 = _helper.GetNodeUV( face, n10.Node(), n11.Node() );
1039   gp_XY uv11 = _helper.GetNodeUV( face, n11.Node(), n10.Node() );
1040   gp_XY uv1  = uv10 * ( 1 - u1 ) + uv11 * u1;
1041
1042   gp_XY uv20 = _helper.GetNodeUV( face, n20.Node(), n21.Node() );
1043   gp_XY uv21 = _helper.GetNodeUV( face, n21.Node(), n20.Node() );
1044   gp_XY uv2  = uv20 * ( 1 - u2 ) + uv21 * u2;
1045
1046   intUV = 0.5 * ( uv1 + uv2 );
1047
1048   // compute precise UV and XYZ by projecting intPnt to the FACE
1049
1050   Handle(ShapeAnalysis_Surface) surface = _helper.GetSurface( face );
1051   intUV = surface->NextValueOfUV( intUV, intPnt, _tol );
1052   if ( surface->Gap() > Min( len1, len2 ))
1053     intUV = surface->ValueOfUV( intPnt, _tol );
1054
1055   intPnt = surface->Value( intUV );
1056
1057   return true;
1058 }
1059
1060 //================================================================================
1061 /*!
1062  * \brief Setup predicates to detect nodes on FACE boundary
1063  *  \param [in] shape - shape containing FACEs to mesh
1064  */
1065 //================================================================================
1066
1067 void BLSURFPlugin_EnforcedMesh1D::setupPredicates( const TopoDS_Shape& theShape )
1068 {
1069   if ( _onVertexPredicate )
1070     return;
1071
1072   _onVertexPredicate.reset( new TPredicate() );
1073   _onVertexPredicate->SetTolerance( _tol );
1074   _onVertexPredicate->SetMesh( _helper.GetMeshDS() );
1075   {
1076     TopTools_IndexedMapOfShape vertices;
1077     TopExp::MapShapes( theShape, TopAbs_VERTEX, vertices );
1078     TopoDS_Compound vCompound;
1079     BRep_Builder    builder;
1080     builder.MakeCompound( vCompound );
1081     for ( const TopoDS_Shape& v : vertices )
1082       builder.Add( vCompound, v );
1083
1084     _onVertexPredicate->SetShape( vCompound, SMDSAbs_Node );
1085   }
1086
1087
1088   _onEdgePredicate.reset( new TPredicate() );
1089   _onEdgePredicate->SetTolerance( _tol );
1090   _onEdgePredicate->SetMesh( _helper.GetMeshDS() );
1091   {
1092     TopTools_IndexedMapOfShape edges;
1093     TopExp::MapShapes( theShape, TopAbs_EDGE, edges );
1094     TopoDS_Compound eCompound;
1095     BRep_Builder    builder;
1096     builder.MakeCompound( eCompound );
1097     for ( const TopoDS_Shape& e : edges )
1098       builder.Add( eCompound, e );
1099
1100     _onEdgePredicate->SetShape( eCompound, SMDSAbs_Node );
1101   }
1102   return;
1103 }
1104
1105 //================================================================================
1106 /*!
1107  * \brief Find or create a group of edges with given name
1108  */
1109 //================================================================================
1110
1111 SMDS_MeshGroup* BLSURFPlugin_EnforcedMesh1D::getGroup( const std::string& groupName )
1112 {
1113   SMDS_MeshGroup* group = nullptr;
1114   if ( !groupName.empty() )
1115   {
1116     if ( _name2Group.count( groupName ))
1117       return _name2Group[ groupName ];
1118
1119     // find existing group
1120     for ( SMESH_Mesh::GroupIteratorPtr grIt = _mesh->GetGroups(); grIt->more(); )
1121     {
1122       SMESH_Group*     grp = grIt->next();
1123       SMESHDS_Group* grpDS = dynamic_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1124       if ( grpDS &&
1125            grpDS->GetType() == SMDSAbs_Edge &&
1126            groupName == grp->GetName() )
1127       {
1128         _name2Group[ groupName ] = & grpDS->SMDSGroup();
1129         return & grpDS->SMDSGroup();
1130       }
1131     }
1132
1133     // create a new group
1134     SMESH_Group*     grp = _mesh->AddGroup( SMDSAbs_Edge, groupName.c_str() );
1135     SMESHDS_Group* grpDS = static_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1136
1137     group = & grpDS->SMDSGroup();
1138     _name2Group[ groupName ] = group;
1139   }
1140   return group;
1141 }
1142
1143 //================================================================================
1144 /*!
1145  * \brief Look for a node dividing a given EDGE
1146  */
1147 //================================================================================
1148
1149 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::findNodeOnEdge( const gp_Pnt&      p,
1150                                                                   const TopoDS_Edge& edge )
1151 {
1152   // look for an equal node on the EDGE
1153   if ( std::vector< const SMDS_MeshNode* >* nodesOnE = _nodesOnEdge.ChangeSeek( edge ))
1154     for ( const SMDS_MeshNode* n : *nodesOnE )
1155       if ( p.SquareDistance( SMESH_NodeXYZ( n )) < _tol * _tol )
1156       {
1157         return n;
1158       }
1159
1160   return nullptr;
1161 }
1162
1163 //================================================================================
1164 /*!
1165  * \brief Add a node to an EDGE
1166  */
1167 //================================================================================
1168
1169 void BLSURFPlugin_EnforcedMesh1D::addNodeOnEdge( const SMDS_MeshNode* node,
1170                                                  const TopoDS_Edge&   edge,
1171                                                  const double         u)
1172 {
1173   _mesh->GetMeshDS()->SetNodeOnEdge( node, edge, u );
1174
1175   std::vector< const SMDS_MeshNode* > * nodesOnE = _nodesOnEdge.ChangeSeek( edge );
1176   if ( !nodesOnE )
1177     nodesOnE = _nodesOnEdge.Bound( edge, std::vector< const SMDS_MeshNode* >() );
1178
1179   nodesOnE->push_back( node );
1180 }
1181
1182 //================================================================================
1183 /*!
1184  * \brief Create EDGEs by dividing a given EDGE by nodes on it
1185  *  \param [in] edge - the EDGE to divide
1186  *  \param [in] nodes - the nodes to divide by
1187  */
1188 //================================================================================
1189
1190 void BLSURFPlugin_EnforcedMesh1D::
1191
1192 splitEdgeByNodes( TopoDS_Edge edge, const std::vector< const SMDS_MeshNode* >& nodes )
1193 {
1194   if ( nodes.empty() )
1195     return;
1196
1197   edge.Orientation( TopAbs_FORWARD );
1198
1199   TopoDS_Vertex v0 = _helper.IthVertex( 0, edge );
1200   TopoDS_Vertex v1 = _helper.IthVertex( 1, edge );
1201
1202   // create VERTEXes and sort them along the EDGE
1203
1204   std::map< double, TopoDS_Vertex > sortedVertices;
1205
1206   BRepAdaptor_Curve curve( edge );
1207   for ( SMESH_NodeXYZ pn : nodes )
1208   {
1209     gp_Pnt projPnt;
1210     double u;
1211     ShapeAnalysis_Curve().Project( curve, pn, _tol, projPnt, u, false );
1212     projPnt = curve.Value( u );
1213
1214     TopoDS_Vertex v = BRepBuilderAPI_MakeVertex( projPnt );
1215
1216     sortedVertices.insert({ u, v });
1217
1218     _nOnE2Vertex[ pn.Node() ] = v;
1219   }
1220   sortedVertices.insert({ BRep_Tool::Parameter( v1, edge ), v1 });
1221
1222
1223   // create EDGEs
1224
1225   BRep_Builder builder;
1226   std::vector< TopoDS_Edge >& newEdges = *_edgeSplitsOfEdge.Bound( edge, TEdge2Edges::value_type());
1227
1228   double u0 = BRep_Tool::Parameter( v0, edge );
1229   for ( auto& u2v : sortedVertices )
1230   {
1231     double u1 = u2v.first;
1232     v1        = u2v.second;
1233
1234     TopoDS_Shape newShape = edge.EmptyCopied();
1235     TopoDS_Edge  newEdge  = TopoDS::Edge( newShape );
1236     builder.Add( newEdge, v0 );
1237     builder.Add( newEdge, v1 );
1238     builder.Range( newEdge, u0, u1 );
1239     newEdges.push_back( newEdge );
1240
1241     v0 = v1;
1242     u0 = u1;
1243   }
1244
1245   return;
1246 }
1247
1248 //================================================================================
1249 /*!
1250  * \brief  Return true if there are enforced segments on a FACE. Start iteration on segments
1251  */
1252 //================================================================================
1253
1254 bool BLSURFPlugin_EnforcedMesh1D::HasSegmentsOnFace( const TopoDS_Face& face )
1255 {
1256   _segIterator.reset();
1257
1258   if ( _faces.Contains( face ))
1259   {
1260     _currentFace = face;
1261     _segIterator = _helper.GetMeshDS()->MeshElements( face )->GetElements();
1262
1263     _helper.SetSubShape( face );
1264
1265     return _segIterator->more();
1266   }
1267   return false;
1268 }
1269
1270
1271 //================================================================================
1272 /*!
1273  * \brief Return next segment on the FACE
1274  */
1275 //================================================================================
1276
1277 bool BLSURFPlugin_EnforcedMesh1D::NextSegment( Segmemnt &                   seg,
1278                                                TopTools_IndexedMapOfShape & vertexTags )
1279 {
1280   if ( _segIterator && _segIterator->more() )
1281   {
1282     const SMDS_MeshElement* segment = _segIterator->next();
1283
1284     while ( segment->GetType() != SMDSAbs_Edge && _segIterator->more() )
1285       segment = _segIterator->next();
1286     if ( segment->GetType() != SMDSAbs_Edge )
1287       return false;
1288
1289     const SMDS_MeshNode * node[2] = { segment->GetNode(0),
1290                                       segment->GetNode(1) };
1291
1292     seg._tag = _segTag++;
1293
1294     seg._xyz[0] = SMESH_NodeXYZ( node[0]);
1295     seg._xyz[1] = SMESH_NodeXYZ( node[1]);
1296
1297     seg._uv[0] = _helper.GetNodeUV( _currentFace, node[0], node[1] );
1298     seg._uv[1] = _helper.GetNodeUV( _currentFace, node[1], node[0] );
1299
1300     seg._pcurve = new Geom2d_Line( seg._uv[0], ( seg._uv[1] - seg._uv[0] ));
1301
1302     seg._u[0] = 0.0;
1303     seg._u[1] = ( seg._uv[1] - seg._uv[0] ).Modulus();
1304
1305     for ( int i = 0; i < 2; ++i )
1306     {
1307       auto n2v = _nOnE2Vertex.find( node[i] ); // find VERTEX by node
1308
1309       if ( n2v != _nOnE2Vertex.end() )
1310         seg._vTag[i] = vertexTags.Add( n2v->second );
1311       else
1312         seg._vTag[i] = _nodeTag0 + _nodeTags.Add( node[i] );
1313     }
1314
1315     if ( !_isQuadratic )
1316       _mesh->GetMeshDS()->UnSetMeshElementOnShape( segment, _currentFace );
1317
1318     return true;
1319   }
1320   return false;
1321 }
1322
1323 //================================================================================
1324 /*!
1325  * \brief Return enforced node by the tag that was returned by Segmemnt::_vTag[i]
1326  */
1327 //================================================================================
1328
1329 const SMDS_MeshNode*
1330 BLSURFPlugin_EnforcedMesh1D::GetNodeByTag( int                                tag,
1331                                            const TopTools_IndexedMapOfShape & vertexTags )
1332 {
1333   const SMDS_MeshNode* node = nullptr;
1334
1335   if ( tag <= vertexTags.Size() && !_nOnE2Vertex.empty() )
1336   {
1337     const TopoDS_Shape& vertex = vertexTags( tag );
1338
1339     for ( auto n2v = _nOnE2Vertex.begin(); n2v != _nOnE2Vertex.end(); ++n2v )
1340       if ( vertex.IsSame( n2v->second ))
1341       {
1342         node = n2v->first;
1343         _nOnE2Vertex.erase( n2v );
1344         return node;
1345       }
1346   }
1347
1348   tag -= _nodeTag0;
1349
1350   bool isValid = ( 0 < tag && tag <= _nodeTags.Size() );
1351   if ( isValid )
1352     node = _nodeTags( tag );
1353
1354   return node;
1355 }
1356
1357 //================================================================================
1358 /*!
1359  * \brief Return true if a tag corresponds to the tag of enforced segment
1360  *  that was returned by Segment::_tag
1361  */
1362 //================================================================================
1363
1364 bool BLSURFPlugin_EnforcedMesh1D::IsSegmentTag( int tag ) const
1365 {
1366   return ( _segTag0 <= tag && tag <= _segTag );
1367 }
1368
1369 //================================================================================
1370 /*!
1371  * \brief Return tag of EDGE by tags of its splits
1372  */
1373 //================================================================================
1374
1375 int BLSURFPlugin_EnforcedMesh1D::GetTagOfSplitEdge( int splitTag ) const
1376 {
1377   auto sTag2eTag = _splitTags2EdgeTag.find( splitTag );
1378   return ( sTag2eTag == _splitTags2EdgeTag.end() ) ? splitTag : sTag2eTag->second;
1379 }