]> SALOME platform Git repositories - plugins/blsurfplugin.git/blob - src/BLSURFPlugin/BLSURFPlugin_EnforcedMesh1D.cxx
Salome HOME
Merge branch 'V9_9_BR'
[plugins/blsurfplugin.git] / src / BLSURFPlugin / BLSURFPlugin_EnforcedMesh1D.cxx
1 // Copyright (C) 2007-2022  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   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 }