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