src/BLSURFPlugin/BLSURFPlugin_EnforcedMesh1D.cxx

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