1 // Copyright (C) 2007-2021 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // File : BLSURFPlugin_EnforcedMesh1D.cxx
20 // Author : Edward AGAPOV (OCC)
22 #include "BLSURFPlugin_EnforcedMesh1D.hxx"
24 #include "BLSURFPlugin_BLSURF.hxx"
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>
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>
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(); }
56 //================================================================================
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
67 //================================================================================
69 const SMDS_MeshNode* findNode( const gp_Pnt& p,
71 const SMDS_MeshElement* segment,
72 BLSURFPlugin_EnforcedMesh1D::TNodesOnSeg & segInternalNodes,
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();
83 auto seg2nodes = segInternalNodes.find( segment );
84 if ( seg2nodes != segInternalNodes.end() )
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 ))
97 //================================================================================
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
107 //================================================================================
109 void orientSegments( const std::vector< const SMDS_MeshElement* > & braSegs,
110 const std::vector< const SMDS_MeshNode* > & braNodes,
111 const size_t nodeIndex,
114 const SMDS_MeshElement* toReverse[2] = { nullptr, nullptr };
116 if ( nodeIndex > 0 &&
117 braSegs[ nodeIndex - 1 ]->GetNode(1) != braNodes[ nodeIndex ])
118 toReverse[ 0 ] = braSegs[ nodeIndex - 1 ];
120 if ( nodeIndex < braSegs.size() &&
121 braSegs[ nodeIndex ]->GetNode(0) != braNodes[ nodeIndex ])
122 toReverse[ bool( toReverse[0]) ] = braSegs[ nodeIndex ];
127 SMESH_MeshEditor editor( mesh );
128 for ( int i = 0; i < 2; ++i )
130 editor.Reorient( toReverse[ i ]);
135 //================================================================================
137 * \brief Create enforced mesh segments in a mesh
138 * \param [in] helper - contains the mesh and the shape to compute
139 * \param [in] hyp - hypothesis containing data of enforced meshes
141 //================================================================================
143 BLSURFPlugin_EnforcedMesh1D::BLSURFPlugin_EnforcedMesh1D( SMESH_MesherHelper& helper,
144 const BLSURFPlugin_Hypothesis* hyp )
145 : _mesh ( helper.GetMesh() ),
146 _shape( helper.GetSubShape() ),
148 _isQuadratic( helper.GetIsQuadratic() )
150 if ( !hyp || !_mesh || hyp->GetEnforcedMeshes().empty() )
153 _tol = 2 * BRep_Tool::MaxTolerance( _shape, TopAbs_VERTEX );
154 _segTag = 1000 + helper.Count( _shape, TopAbs_EDGE, /*ignoreSame=*/false );
156 _nodeTag0 = 1000 + helper.Count( _shape, TopAbs_VERTEX, /*ignoreSame=*/false );
158 for ( const BLSURFPlugin_Hypothesis::EnforcedMesh& enfMesh : hyp->GetEnforcedMeshes() )
160 copyEnforcedMesh( enfMesh, hyp, _shape );
163 for ( const SMDS_MeshElement* segment : _segmentsOnSeveralShapes )
165 splitSegmentOnSeveralShapes( segment );
168 for ( const TopoDS_Shape & face : _faces )
170 splitSelfIntersectingSegments( face );
173 for ( TEdge2Nodes::Iterator e2nn( _nodesOnEdge ); e2nn.More(); e2nn.Next() )
175 splitEdgeByNodes( e2nn.Key(), e2nn.Value() );
179 //================================================================================
181 * \brief Convert enforced segments to quadratic
183 //================================================================================
185 BLSURFPlugin_EnforcedMesh1D::~BLSURFPlugin_EnforcedMesh1D()
190 _helper.SetIsQuadratic( true );
192 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
193 std::vector<const SMDS_MeshNode *> nodes(2);
194 std::vector<const SMDS_MeshElement *> trias;
196 SMDS_MeshGroup* group = nullptr;
197 if ( !_name2Group.empty() )
198 group = _name2Group.begin()->second;
200 for ( const TopoDS_Shape& f : _faces )
201 if ( HasSegmentsOnFace( TopoDS::Face( f )))
202 while ( _segIterator->more() )
204 const SMDS_MeshElement* segment = _segIterator->next();
205 if ( segment->GetType() != SMDSAbs_Edge )
208 nodes[0] = segment->GetNode(0);
209 nodes[1] = segment->GetNode(1);
210 if ( meshDS->GetElementsByNodes( nodes, trias, SMDSAbs_Face ))
212 if ( !group || !group->Contains( segment ))
214 SMDS_MeshGroup* otherGroup = nullptr;
215 if ( _name2Group.size() > bool( group ))
216 for ( auto & n2g : _name2Group )
217 if ( n2g.second != group && n2g.second->Contains( segment ))
219 otherGroup = n2g.second;
225 _helper.AddTLinks( static_cast< const SMDS_MeshFace* >( trias[0] ));
227 meshDS->RemoveFreeElement( segment, meshDS->MeshElements( f ), /*fromGroup=*/false );
229 SMDS_MeshElement* quadSegment = _helper.AddEdge( nodes[0], nodes[1],
230 /*id=*/0, /*force3d=*/false );
231 if ( group && segment != quadSegment )
233 group->Remove( segment );
234 group->Add( quadSegment );
241 //================================================================================
243 * \brief Return EDGEs resulted from division of FACE boundary by enforced segments
244 * \param [in] edge - possibly divided EDGE
245 * \param [out] splits - split EDGEs
246 * \return bool - true if the EDGE is split
248 //================================================================================
250 bool BLSURFPlugin_EnforcedMesh1D::GetSplitsOfEdge( const TopoDS_Edge& edge,
251 std::vector< TopoDS_Edge > & splits,
252 TopTools_IndexedMapOfShape & edgeTags )
254 std::vector< TopoDS_Edge > * splitsInMap = _edgeSplitsOfEdge.ChangeSeek( edge );
258 splits.insert( splits.end(), splitsInMap->begin(), splitsInMap->end() );
260 int eTag = edgeTags.Add( edge );
262 size_t index = splits.size() - 1;
263 for ( size_t i = 0; i < splitsInMap->size(); ++i, --index )
265 int splitTag = edgeTags.Add( splits[ index ]);
266 _splitTags2EdgeTag[ splitTag ] = eTag;
268 if ( edge.Orientation() == TopAbs_REVERSED )
269 splits[ index ].Reverse();
276 //================================================================================
278 * \brief Add 1D elements to the mesh
280 //================================================================================
282 void BLSURFPlugin_EnforcedMesh1D::
284 copyEnforcedMesh( const BLSURFPlugin_Hypothesis::EnforcedMesh& theEnfMesh,
285 const BLSURFPlugin_Hypothesis* theHyp,
286 const TopoDS_Shape& theShape)
289 SMDS_ElemIteratorPtr segIt = theHyp->GetEnforcedSegments( theEnfMesh, mesh1D );
290 if ( !segIt->more() )
293 // setup predicates to detect nodes on FACE boundary
294 setupPredicates( theShape );
296 SMDS_MeshGroup* group = getGroup( theEnfMesh._groupName );
297 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
299 // get ordered nodes and segments of theEnfMesh
300 SMESH_MeshAlgos::TElemGroupVector edgeBranches;
301 SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
302 SMESH_MeshAlgos::Get1DBranches( segIt, edgeBranches, nodeBranches );
305 // Copy nodes and segments from an enforced mesh to my mesh
307 TopoDS_Shape vertex, edge;
309 // first treat ends of branches that can be shared by branches
310 for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
312 std::vector< const SMDS_MeshNode* > & braNodes = nodeBranches[ iB ];
313 std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
315 for ( int isLast = 0; isLast < 2; ++isLast )
317 const SMDS_MeshNode* enfNode = isLast ? braNodes.back() : braNodes[0];
318 if ( meshDS->Contains( enfNode ))
319 continue; // already in my mesh
321 const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
323 orientSegments( braSegs, braNodes, isLast ? 0 : braNodes.size() - 1, mesh1D );
325 // replace enfNode at branch ends by newNode
326 SMESH_NodeXYZ enfPnt( newNode ? newNode : enfNode );
327 for ( std::vector< const SMDS_MeshNode* > & braNodes : nodeBranches )
329 for ( int isLast = 0; isLast < 2; ++isLast )
331 const SMDS_MeshNode* & endNode = isLast ? braNodes.back() : braNodes[0];
332 if ( endNode == enfNode || enfPnt.SquareDistance( endNode ) < _tol*_tol )
337 } // loop on branch ends
338 } // loop on nodeBranches
340 // copy nodes and segments
342 for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
344 std::vector< const SMDS_MeshNode* > & braNodes = nodeBranches[ iB ];
345 std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
347 // copy nodes of the branch
348 for ( size_t i = 0; i < braNodes.size(); ++i )
350 const SMDS_MeshNode* & enfNode = braNodes[ i ];
351 const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
353 if ( !newNode ) // orient segments to be able to get enforced not projected node
354 orientSegments( braSegs, braNodes, i, mesh1D );
359 // copy segments of the branch
360 for ( size_t i = 0; i < braSegs.size(); ++i )
362 //braSegs[ i ] = nullptr;
364 const SMDS_MeshNode* node0 = braNodes[ i ];
365 const SMDS_MeshNode* node1 = braNodes[ i + 1 ];
366 if ( !node0 && !node1 )
369 TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
370 TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
371 if ( shape0.IsNull() && shape1.IsNull() )
374 if ( !node0 && shape1.ShapeType() != TopAbs_FACE )
376 if ( !node1 && shape0.ShapeType() != TopAbs_FACE )
379 if ( !node0 || !node1 ) // create missing node at location of enforced node projected nowhere
381 SMESH_NodeXYZ pn = braSegs[i]->GetNode( !node1 );
382 ( node0 ? node1 : node0 ) = _helper.AddNode( pn->X(), pn->Y(), pn->Z() );
385 SMDS_MeshEdge* newSeg = _helper.AddEdge( node0, node1 );
387 group->Add( newSeg );
388 braSegs[ i ] = newSeg;
390 // check if the both nodes are on the same FACE
391 TopoDS_Shape face = shape0;
392 if ( !shape0.IsSame( shape1 ) && !shape0.IsNull() && !shape1.IsNull() )
394 if ( shape0.ShapeType() == TopAbs_FACE &&
395 _helper.IsSubShape( shape1, shape0 ))
399 else if ( shape1.ShapeType() == TopAbs_FACE &&
400 _helper.IsSubShape( shape0, shape1 ))
404 else // try to find a FACE by projecting a segment middle point
407 gp_Pnt middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
408 //BLSURFPlugin_BLSURF::projectionPoint projPnt =
409 BLSURFPlugin_BLSURF::getProjectionPoint( TopoDS::Face( face ), middlePnt );
411 if ( !face.IsNull() &&
412 ( !_helper.IsSubShape( shape0, face ) ||
413 !_helper.IsSubShape( shape1, face ) ))
417 if ( !face.IsNull() && face.ShapeType() == TopAbs_FACE )
419 meshDS->SetMeshElementOnShape( newSeg, face );
423 if ( face.IsNull() || shape0.IsNull() || shape1.IsNull() )
425 _segmentsOnSeveralShapes.push_back( newSeg );
428 } // loop on branch segments
430 } // loop on branches
435 //================================================================================
437 * \brief Create a copy of a node of enforced mesh in my mesh
438 * \param [in] enfNode - enforced node
439 * \return const SMDS_MeshNode* - a node in my mesh
441 //================================================================================
443 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::copyEnforcedNode( const SMDS_MeshNode* enfNode )
445 SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
447 if ( !enfNode || meshDS->Contains( enfNode ))
448 return enfNode; // already in my mesh
450 SMESH_NodeXYZ enfPnt = enfNode;
452 const SMDS_MeshNode* newNode = nullptr;
454 // check if enfNode is on VERTEX
455 TopoDS_Vertex vertex;
456 if ( _onVertexPredicate->IsSatisfy( enfNode, &vertex ))
458 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
459 newNode = SMESH_Algo::VertexNode( vertex, meshDS );
462 // check if enfNode is on EDGE
463 bool setNodeOnEdge = false;
467 setNodeOnEdge = _onEdgePredicate->IsSatisfy( enfNode, &edge );
470 newNode = findNodeOnEdge( enfPnt, edge );
471 setNodeOnEdge = !newNode;
475 // create a new node and set it on FACE
479 BLSURFPlugin_BLSURF::projectionPoint projPnt =
480 BLSURFPlugin_BLSURF::getProjectionPoint( face, enfPnt, /*allowStateON=*/true );
481 if ( face.IsNull() ) return newNode;
483 if ( projPnt.state == TopAbs_ON )
485 SMDS_MeshNode* projNode = const_cast< SMDS_MeshNode* >( enfNode );
486 projNode->setXYZ( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
489 if ( _onVertexPredicate->IsSatisfy( projNode, &vertex ))
491 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
492 newNode = SMESH_Algo::VertexNode( vertex, meshDS );
494 else if (( setNodeOnEdge = _onEdgePredicate->IsSatisfy( projNode, &edge )))
496 newNode = findNodeOnEdge( projPnt.xyz, edge );
497 setNodeOnEdge = !newNode;
499 projNode->setXYZ( enfPnt.X(), enfPnt.Y(), enfPnt.Z() );
503 newNode = meshDS->AddNode( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
505 if ( vertex.IsNull() && edge.IsNull() )
506 meshDS->SetNodeOnFace( newNode, face, projPnt.uv.X(), projPnt.uv.Y() );
509 // set the new node on EDGE
510 if ( newNode && setNodeOnEdge )
512 addNodeOnEdge( newNode, edge );
518 //================================================================================
520 * \brief Split a segment whose nodes are on different FACEs into smaller parts
521 * lying each on one FACE
523 //================================================================================
525 void BLSURFPlugin_EnforcedMesh1D::splitSegmentOnSeveralShapes( const SMDS_MeshElement* segment )
527 const SMDS_MeshNode* node0 = segment->GetNode(0);
528 const SMDS_MeshNode* node1 = segment->GetNode(1);
529 SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
531 TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
532 TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
534 if ( shape0.IsNull() )
536 std::swap( node0, node1 );
541 gp_XYZ xyz0 = SMESH_NodeXYZ( node0 );
542 gp_XYZ xyz1 = SMESH_NodeXYZ( node1 );
543 gp_XYZ segDir = ( xyz1 - xyz0 ).Normalized();
545 //std::map< double, const SMDS_MeshNode* > mediumNodes; // nodes splitting the segment
547 while ( !shape0.IsSame( shape1 )) // move along the segment till shape1
549 if ( shape0.ShapeType() == TopAbs_FACE ) // make a node on an EDGE of the FACE
550 { // ----------------------------------
553 gp_Pnt edgeIntPnt = getEdgeIntersection( shape0, xyz0, xyz1, edge, paramOnE );
557 // check if edgeIntPnt in on VERTEX
558 TopoDS_Vertex vertex;
559 for ( int iV = 0; iV < 2 && vertex.IsNull(); ++iV )
561 TopoDS_Vertex v = _helper.IthVertex( iV, edge );
562 if ( edgeIntPnt.SquareDistance( BRep_Tool::Pnt( v )) < _tol *_tol )
566 // make a node on the EDGE
567 const SMDS_MeshNode* nodeOnEdge = nullptr;
568 if ( vertex.IsNull() )
570 nodeOnEdge = findNodeOnEdge( edgeIntPnt, edge );
573 if ( shape1.IsNull() && node1 )
577 meshDS->MoveNode( nodeOnEdge, edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
581 nodeOnEdge = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
583 addNodeOnEdge( nodeOnEdge, edge, paramOnE );
588 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
589 nodeOnEdge = SMESH_Algo::VertexNode( vertex, meshDS );
592 // create a sub-segment
593 SMDS_MeshElement* newSeg = _helper.AddEdge( node0, nodeOnEdge );
594 meshDS->SetMeshElementOnShape( newSeg, shape0 );
596 SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
599 xyz0 = SMESH_NodeXYZ( node0 );
600 if ( vertex.IsNull() )
606 else // shape0 is EDGE or VERTEX; look for the next FACE
607 { // ------------------------------------------------
609 if ( !shape1.IsNull() &&
610 shape1.ShapeType() == TopAbs_FACE &&
611 _helper.IsSubShape( shape0, shape1 )) // shape0 belongs to FACE shape1
613 SMDS_MeshElement* newSeg = _helper.AddEdge( node0, node1 );
614 SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
615 meshDS->SetMeshElementOnShape( newSeg, shape1 );
620 double shift = 10 * _tol;
621 for ( int nbAttemp = 0; face.IsNull() && nbAttemp < 10; ++nbAttemp )
623 xyz0 += segDir * shift;
625 BLSURFPlugin_BLSURF::getProjectionPoint( face, xyz0 );
627 if ( !face.IsNull() )
629 if ( _helper.IsSubShape( shape0, face ))
641 } //while ( !shape0.IsSame( shape1 ))
643 meshDS->RemoveFreeElement( segment, /*submesh=*/nullptr );
646 //================================================================================
648 * \brief Find intersection of FACE EDGEs and a segment
649 * \param [in] theFace - the FACE
650 * \param [in] thePnt0 - first end of the segment
651 * \param [in] thePnt1 - last end of the segment
652 * \param [out] theFounfEdge - return the intersected EDGE
653 * \param [out] theParamOnEdge - return parameter of intersection point on EDGE
654 * \return gp_XYZ - point on an EDGE closest to the segment
656 //================================================================================
658 gp_Pnt BLSURFPlugin_EnforcedMesh1D::getEdgeIntersection( const TopoDS_Shape& theFaceOrEdge,
659 const gp_XYZ& thePnt0,
660 const gp_XYZ& thePnt1,
661 TopoDS_Edge & theFounfEdge,
662 double & theParamOnEdge)
664 const double segLen = ( thePnt1 - thePnt0 ).Modulus();
665 const double maxSegDist2 = segLen * segLen * 0.5 * 0.5;
667 Handle(Geom_Line) segLine = new Geom_Line( thePnt0, thePnt1 - thePnt0 );
668 GeomAdaptor_Curve segLineAdpt( segLine, 0, segLen );
670 TopTools_MapOfShape edges;
671 double minParamOnSeg = segLen;
673 for ( TopExp_Explorer edgeExp( theFaceOrEdge, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
675 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
676 if ( !edges.Add( edge ))
679 Extrema_ExtCC extrema( segLineAdpt, BRepAdaptor_Curve( edge ), _tol, _tol );
681 if ( extrema.IsDone() && !extrema.IsParallel() )
682 for ( int i = 1, nb = extrema.NbExt(); i <= nb; ++i )
683 if ( extrema.SquareDistance( i ) < maxSegDist2 )
685 Extrema_POnCurv pOnSeg, pOnEdge;
686 extrema.Points( i, pOnSeg, pOnEdge );
687 double paramOnSeg = pOnSeg.Parameter();
688 if ( 0 < paramOnSeg && paramOnSeg < minParamOnSeg )
690 minParamOnSeg = paramOnSeg;
691 foundPnt = pOnEdge.Value();
693 theParamOnEdge = pOnEdge.Parameter();
700 //================================================================================
702 * \brief Split self-intersecting segments on a given FACE
704 //================================================================================
706 void BLSURFPlugin_EnforcedMesh1D::splitSelfIntersectingSegments( const TopoDS_Shape & theFace )
708 const TopoDS_Face& face = TopoDS::Face( theFace );
710 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
711 SMESHDS_SubMesh* sm = meshDS->MeshElements( face );
712 if ( !sm || sm->NbElements() <= 1 )
715 // get ordered nodes and segments on the face
716 SMESH_MeshAlgos::TElemGroupVector edgeBranches;
717 SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
718 SMESH_MeshAlgos::Get1DBranches( sm->GetElements(), edgeBranches, nodeBranches );
720 // create element searcher
721 SMESH_ElementSearcher* elemSearcher;
722 SMESHUtils::Deleter< SMESH_ElementSearcher > elSearchdeleter;
723 std::vector< const SMDS_MeshElement* > foundElems;
725 std::vector< SMDS_ElemIteratorPtr > elemItVec;
726 for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
728 SMDS_ElemIteratorPtr segIt =
729 boost::make_shared< SMDS_ElementVectorIterator >( braSegs.begin(), braSegs.end() );
730 elemItVec.push_back( segIt );
732 typedef SMDS_IteratorOnIterators< const SMDS_MeshElement*,
733 std::vector< SMDS_ElemIteratorPtr > > TVecIterator;
734 SMDS_ElemIteratorPtr segIt = boost::make_shared< TVecIterator >( elemItVec );
736 elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, segIt, _tol );
737 elSearchdeleter._obj = elemSearcher;
739 // force usage of iterators before they die
740 elemSearcher->FindElementsByPoint( gp_Pnt( 0,0,1e+20), SMDSAbs_Edge, foundElems );
744 // Find intersecting segments
746 std::map< const SMDS_MeshElement* , std::vector< const SMDS_MeshNode* > > segInternalNodes;
747 SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
749 for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
751 braSegs.push_back( braSegs.back() );
752 const SMDS_MeshElement* prevSeg = nullptr;
754 for ( size_t i = 0, nb = braSegs.size() - 1; i < nb; ++i )
756 const SMDS_MeshElement* seg = braSegs[ i ];
757 const SMDS_MeshElement* nextSeg = braSegs[ i + 1 ];
758 const SMDS_MeshNode* node0 = seg->GetNode(0);
759 const SMDS_MeshNode* node1 = seg->GetNode(1);
761 gp_XYZ xyz0 = SMESH_NodeXYZ( node0 );
762 gp_XYZ xyz1 = SMESH_NodeXYZ( node1 );
763 gp_XYZ middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
764 double segLen = ( xyz0 - xyz1 ).Modulus();
767 elemSearcher->GetElementsInSphere( middlePnt, 0.5 * segLen + _tol, SMDSAbs_Edge, foundElems );
769 for ( const SMDS_MeshElement* closeSeg : foundElems )
771 if ( closeSeg == prevSeg ||
773 closeSeg == nextSeg )
778 if ( intersectSegments( seg, closeSeg, face, intPnt, uv ))
781 const SMDS_MeshNode* intNode0 = findNode( intPnt, _tol, seg, segInternalNodes, i0 );
782 const SMDS_MeshNode* intNode1 = findNode( intPnt, _tol, closeSeg, segInternalNodes, i1 );
783 if ( !intNode0 && intNode1 )
785 segInternalNodes[ seg ].push_back( intNode1 );
787 else if ( !intNode1 && intNode0 )
789 segInternalNodes[ closeSeg ].push_back( intNode0 );
791 else if ( intNode1 && intNode0 )
793 if ( intNode1 == intNode0 )
795 if ( i0 < 0 && i1 < 0 )
797 nodeGroupsToMerge.push_back // merge end nodes
798 ( std::list< const SMDS_MeshNode* >({ intNode0, intNode1 }));
802 segInternalNodes[ closeSeg ][ i1 ] = intNode0;
806 segInternalNodes[ seg ][ i0 ] = intNode1;
808 else // two internal nodes coincide
810 segInternalNodes[ seg ][ i0 ] = intNode1;
811 nodeGroupsToMerge.push_back
812 ( std::list< const SMDS_MeshNode* >({ intNode1, intNode0 }));
815 else // ( !intNode1 && !intNode0 )
817 intNode0 = _helper.AddNode( intPnt.X(), intPnt.Y(), intPnt.Z() );
818 meshDS->SetNodeOnFace( intNode0, face, uv.X(), uv.Y() );
819 segInternalNodes[ seg ].push_back( intNode0 );
820 segInternalNodes[ closeSeg ].push_back( intNode0 );
827 } // loop on segments of a branch
828 } // loop on branches
831 findIntersectionWithSeamEdge( face, segInternalNodes ); // on periodic FACE
836 for ( auto& seg2nodes : segInternalNodes )
838 const SMDS_MeshElement* seg = seg2nodes.first;
839 std::vector< const SMDS_MeshNode* > & nodes = seg2nodes.second;
840 if ( nodes.empty() ) continue;
842 const SMDS_MeshNode* n0 = seg->GetNode( 0 );
843 const SMDS_MeshNode* n1 = seg->GetNode( 1 );
844 nodes.push_back( n1 );
846 // sort nodes on the segment
847 gp_Pnt p0 = SMESH_NodeXYZ( n0 );
848 std::map< double, const SMDS_MeshNode* > sortedNodes;
849 for ( SMESH_NodeXYZ pn : nodes )
850 sortedNodes.insert({ p0.SquareDistance( pn ), pn.Node() });
853 for ( auto & d2n : sortedNodes )
856 SMDS_MeshElement* newSeg = _helper.AddEdge( n0, n1 );
858 meshDS->SetMeshElementOnShape( newSeg, face );
859 SMESH_MeshEditor::AddToSameGroups( newSeg, seg, meshDS );
861 meshDS->RemoveFreeElement( seg, /*submesh=*/nullptr );
866 SMESH_MeshEditor( _mesh ).MergeNodes( nodeGroupsToMerge );
869 //================================================================================
871 * \brief Find intersections of segments with a seam EDGE on a periodic FACE
873 //================================================================================
875 void BLSURFPlugin_EnforcedMesh1D::findIntersectionWithSeamEdge( const TopoDS_Face & face,
876 TNodesOnSeg & segInternalNodes )
878 _helper.SetSubShape( face );
879 if ( !_helper.HasSeam() )
882 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
883 SMESHDS_SubMesh* sm = meshDS->MeshElements( face );
884 if ( !sm || sm->NbElements() == 0 )
887 TopTools_MapOfShape treatedEdges;
888 for ( TopExp_Explorer edgeExp( face, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
890 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
891 if ( !_helper.IsSeamShape( edge ) || !treatedEdges.Add( edge ))
894 for ( SMDS_ElemIteratorPtr setIt = sm->GetElements(); setIt->more(); )
896 const SMDS_MeshElement* segment = setIt->next();
897 const SMESH_NodeXYZ pn0 = segment->GetNode( 0 );
898 const SMESH_NodeXYZ pn1 = segment->GetNode( 1 );
900 gp_XY uv0 = _helper.GetNodeUV( face, pn0.Node() );
901 gp_XY uv1 = _helper.GetNodeUV( face, pn1.Node() );
903 for ( int iCoo = 1; iCoo <= 2; ++iCoo )
904 if ( iCoo & _helper.GetPeriodicIndex() )
906 double distParam = Abs( uv0.Coord( iCoo ) - uv1.Coord( iCoo ));
907 if ( distParam > 0.5 * _helper.GetPeriod( iCoo ))
909 double paramOnE; TopoDS_Edge edge2;
910 gp_Pnt edgeIntPnt = getEdgeIntersection( edge, pn0, pn1, edge2, paramOnE );
911 if ( edge2.IsNull() )
914 const SMDS_MeshNode* nodeOnSeam = findNodeOnEdge( edgeIntPnt, edge );
915 bool isOnEdge = nodeOnSeam;
917 int indexOnSegment = 3;
919 nodeOnSeam = findNode( edgeIntPnt, _tol, segment, segInternalNodes, indexOnSegment );
923 nodeOnSeam = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
928 addNodeOnEdge( nodeOnSeam, edge, paramOnE );
930 if ( indexOnSegment == 3 )
931 segInternalNodes[ segment ].push_back( nodeOnSeam );
938 //================================================================================
940 * \brief Intersect two segments
941 * \param [in] seg1 - segment 1
942 * \param [in] seg2 - segment 2
943 * \param [in] face - the FACE on which segments lie
944 * \param [out] intPnt - intersection point
945 * \param [out] intUV - UV of intersection point on the FACE
946 * \return bool - true if intersection found
948 //================================================================================
950 bool BLSURFPlugin_EnforcedMesh1D::intersectSegments( const SMDS_MeshElement* seg1,
951 const SMDS_MeshElement* seg2,
952 const TopoDS_Face& face,
954 gp_Pnt2d & intUV ) const
956 SMESH_NodeXYZ n10 = seg1->GetNode(0);
957 SMESH_NodeXYZ n11 = seg1->GetNode(1);
958 SMESH_NodeXYZ n20 = seg2->GetNode(0);
959 SMESH_NodeXYZ n21 = seg2->GetNode(1);
960 if ( n10 == n20 || n10 == n21 || n11 == n20 || n10 == n21 )
963 gp_Lin lin1( n10, n11 - n10 );
964 gp_Lin lin2( n20, n21 - n20 );
966 Extrema_ExtElC extrema( lin1, lin2, Precision::Angular() );
968 if ( !extrema.IsDone() || extrema.IsParallel() )
971 Extrema_POnCurv poc1, poc2;
972 extrema.Points( 1, poc1, poc2 );
974 double len1 = lin1.Direction().XYZ() * ( n11 - n10 );
975 double len2 = lin2.Direction().XYZ() * ( n21 - n20 );
976 double u1 = poc1.Parameter();
977 double u2 = poc2.Parameter();
979 if ( u1 < -_tol || u1 > len1 + _tol )
981 if ( u2 < -_tol || u2 > len2 + _tol )
984 intPnt = 0.5 * ( poc1.Value().XYZ() + poc2.Value().XYZ() );
986 // compute approximate UV
991 gp_XY uv10 = _helper.GetNodeUV( face, n10.Node(), n11.Node() );
992 gp_XY uv11 = _helper.GetNodeUV( face, n11.Node(), n10.Node() );
993 gp_XY uv1 = uv10 * ( 1 - u1 ) + uv11 * u1;
995 gp_XY uv20 = _helper.GetNodeUV( face, n20.Node(), n21.Node() );
996 gp_XY uv21 = _helper.GetNodeUV( face, n21.Node(), n20.Node() );
997 gp_XY uv2 = uv20 * ( 1 - u2 ) + uv21 * u2;
999 intUV = 0.5 * ( uv1 + uv2 );
1001 // compute precise UV and XYZ by projecting intPnt to the FACE
1003 Handle(ShapeAnalysis_Surface) surface = _helper.GetSurface( face );
1004 intUV = surface->NextValueOfUV( intUV, intPnt, _tol );
1005 if ( surface->Gap() > Min( len1, len2 ))
1006 intUV = surface->ValueOfUV( intPnt, _tol );
1008 intPnt = surface->Value( intUV );
1013 //================================================================================
1015 * \brief Setup predicates to detect nodes on FACE boundary
1016 * \param [in] shape - shape containing FACEs to mesh
1018 //================================================================================
1020 void BLSURFPlugin_EnforcedMesh1D::setupPredicates( const TopoDS_Shape& theShape )
1022 if ( _onVertexPredicate )
1025 _onVertexPredicate.reset( new TPredicate() );
1026 _onVertexPredicate->SetTolerance( _tol );
1027 _onVertexPredicate->SetMesh( _helper.GetMeshDS() );
1029 TopTools_IndexedMapOfShape vertices;
1030 TopExp::MapShapes( theShape, TopAbs_VERTEX, vertices );
1031 TopoDS_Compound vCompound;
1032 BRep_Builder builder;
1033 builder.MakeCompound( vCompound );
1034 for ( const TopoDS_Shape& v : vertices )
1035 builder.Add( vCompound, v );
1037 _onVertexPredicate->SetShape( vCompound, SMDSAbs_Node );
1041 _onEdgePredicate.reset( new TPredicate() );
1042 _onEdgePredicate->SetTolerance( _tol );
1043 _onEdgePredicate->SetMesh( _helper.GetMeshDS() );
1045 TopTools_IndexedMapOfShape edges;
1046 TopExp::MapShapes( theShape, TopAbs_EDGE, edges );
1047 TopoDS_Compound eCompound;
1048 BRep_Builder builder;
1049 builder.MakeCompound( eCompound );
1050 for ( const TopoDS_Shape& e : edges )
1051 builder.Add( eCompound, e );
1053 _onEdgePredicate->SetShape( eCompound, SMDSAbs_Node );
1058 //================================================================================
1060 * \brief Find or create a group of edges with given name
1062 //================================================================================
1064 SMDS_MeshGroup* BLSURFPlugin_EnforcedMesh1D::getGroup( const std::string& groupName )
1066 SMDS_MeshGroup* group = nullptr;
1067 if ( !groupName.empty() )
1069 if ( _name2Group.count( groupName ))
1070 return _name2Group[ groupName ];
1072 // find existing group
1073 for ( SMESH_Mesh::GroupIteratorPtr grIt = _mesh->GetGroups(); grIt->more(); )
1075 SMESH_Group* grp = grIt->next();
1076 SMESHDS_Group* grpDS = dynamic_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1078 grpDS->GetType() == SMDSAbs_Edge &&
1079 groupName == grp->GetName() )
1081 _name2Group[ groupName ] = & grpDS->SMDSGroup();
1082 return & grpDS->SMDSGroup();
1086 // create a new group
1087 SMESH_Group* grp = _mesh->AddGroup( SMDSAbs_Edge, groupName.c_str() );
1088 SMESHDS_Group* grpDS = static_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1090 group = & grpDS->SMDSGroup();
1091 _name2Group[ groupName ] = group;
1096 //================================================================================
1098 * \brief Look for a node dividing a given EDGE
1100 //================================================================================
1102 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::findNodeOnEdge( const gp_Pnt& p,
1103 const TopoDS_Edge& edge )
1105 // look for an equal node on the EDGE
1106 if ( std::vector< const SMDS_MeshNode* >* nodesOnE = _nodesOnEdge.ChangeSeek( edge ))
1107 for ( const SMDS_MeshNode* n : *nodesOnE )
1108 if ( p.SquareDistance( SMESH_NodeXYZ( n )) < _tol * _tol )
1116 //================================================================================
1118 * \brief Add a node to an EDGE
1120 //================================================================================
1122 void BLSURFPlugin_EnforcedMesh1D::addNodeOnEdge( const SMDS_MeshNode* node,
1123 const TopoDS_Edge& edge,
1126 _mesh->GetMeshDS()->SetNodeOnEdge( node, edge, u );
1128 std::vector< const SMDS_MeshNode* > * nodesOnE = _nodesOnEdge.ChangeSeek( edge );
1130 nodesOnE = _nodesOnEdge.Bound( edge, std::vector< const SMDS_MeshNode* >() );
1132 nodesOnE->push_back( node );
1135 //================================================================================
1137 * \brief Create EDGEs by dividing a given EDGE by nodes on it
1138 * \param [in] edge - the EDGE to divide
1139 * \param [in] nodes - the nodes to divide by
1141 //================================================================================
1143 void BLSURFPlugin_EnforcedMesh1D::
1145 splitEdgeByNodes( TopoDS_Edge edge, const std::vector< const SMDS_MeshNode* >& nodes )
1147 if ( nodes.empty() )
1150 edge.Orientation( TopAbs_FORWARD );
1152 TopoDS_Vertex v0 = _helper.IthVertex( 0, edge );
1153 TopoDS_Vertex v1 = _helper.IthVertex( 1, edge );
1155 // create VERTEXes and sort them along the EDGE
1157 std::map< double, TopoDS_Vertex > sortedVertices;
1159 BRepAdaptor_Curve curve( edge );
1160 for ( SMESH_NodeXYZ pn : nodes )
1164 ShapeAnalysis_Curve().Project( curve, pn, _tol, projPnt, u, false );
1165 projPnt = curve.Value( u );
1167 TopoDS_Vertex v = BRepBuilderAPI_MakeVertex( projPnt );
1169 sortedVertices.insert({ u, v });
1171 _nOnE2Vertex[ pn.Node() ] = v;
1173 sortedVertices.insert({ BRep_Tool::Parameter( v1, edge ), v1 });
1178 BRep_Builder builder;
1179 std::vector< TopoDS_Edge >& newEdges = *_edgeSplitsOfEdge.Bound( edge, TEdge2Edges::value_type());
1181 double u0 = BRep_Tool::Parameter( v0, edge );
1182 for ( auto& u2v : sortedVertices )
1184 double u1 = u2v.first;
1187 TopoDS_Shape newShape = edge.EmptyCopied();
1188 TopoDS_Edge newEdge = TopoDS::Edge( newShape );
1189 builder.Add( newEdge, v0 );
1190 builder.Add( newEdge, v1 );
1191 builder.Range( newEdge, u0, u1 );
1192 newEdges.push_back( newEdge );
1201 //================================================================================
1203 * \brief Return true if there are enforced segments on a FACE. Start iteration on segments
1205 //================================================================================
1207 bool BLSURFPlugin_EnforcedMesh1D::HasSegmentsOnFace( const TopoDS_Face& face )
1209 _segIterator.reset();
1211 if ( _faces.Contains( face ))
1213 _currentFace = face;
1214 _segIterator = _helper.GetMeshDS()->MeshElements( face )->GetElements();
1216 _helper.SetSubShape( face );
1218 return _segIterator->more();
1224 //================================================================================
1226 * \brief Return next segment on the FACE
1228 //================================================================================
1230 bool BLSURFPlugin_EnforcedMesh1D::NextSegment( Segmemnt & seg,
1231 TopTools_IndexedMapOfShape & vertexTags )
1233 if ( _segIterator && _segIterator->more() )
1235 const SMDS_MeshElement* segment = _segIterator->next();
1237 while ( segment->GetType() != SMDSAbs_Edge && _segIterator->more() )
1238 segment = _segIterator->next();
1239 if ( segment->GetType() != SMDSAbs_Edge )
1242 const SMDS_MeshNode * node[2] = { segment->GetNode(0),
1243 segment->GetNode(1) };
1245 seg._tag = _segTag++;
1247 seg._xyz[0] = SMESH_NodeXYZ( node[0]);
1248 seg._xyz[1] = SMESH_NodeXYZ( node[1]);
1250 seg._uv[0] = _helper.GetNodeUV( _currentFace, node[0], node[1] );
1251 seg._uv[1] = _helper.GetNodeUV( _currentFace, node[1], node[0] );
1253 seg._pcurve = new Geom2d_Line( seg._uv[0], ( seg._uv[1] - seg._uv[0] ));
1256 seg._u[1] = ( seg._uv[1] - seg._uv[0] ).Modulus();
1258 for ( int i = 0; i < 2; ++i )
1260 auto n2v = _nOnE2Vertex.find( node[i] ); // find VERTEX by node
1262 if ( n2v != _nOnE2Vertex.end() )
1263 seg._vTag[i] = vertexTags.Add( n2v->second );
1265 seg._vTag[i] = _nodeTag0 + _nodeTags.Add( node[i] );
1268 if ( !_isQuadratic )
1269 _mesh->GetMeshDS()->UnSetMeshElementOnShape( segment, _currentFace );
1276 //================================================================================
1278 * \brief Return enforced node by the tag that was returned by Segmemnt::_vTag[i]
1280 //================================================================================
1282 const SMDS_MeshNode*
1283 BLSURFPlugin_EnforcedMesh1D::GetNodeByTag( int tag,
1284 const TopTools_IndexedMapOfShape & vertexTags )
1286 const SMDS_MeshNode* node = nullptr;
1288 if ( tag <= vertexTags.Size() && !_nOnE2Vertex.empty() )
1290 const TopoDS_Shape& vertex = vertexTags( tag );
1292 for ( auto n2v = _nOnE2Vertex.begin(); n2v != _nOnE2Vertex.end(); ++n2v )
1293 if ( vertex.IsSame( n2v->second ))
1296 _nOnE2Vertex.erase( n2v );
1303 bool isValid = ( 0 < tag && tag <= _nodeTags.Size() );
1305 node = _nodeTags( tag );
1310 //================================================================================
1312 * \brief Return true if a tag corresponds to the tag of enforced segment
1313 * that was returned by Segment::_tag
1315 //================================================================================
1317 bool BLSURFPlugin_EnforcedMesh1D::IsSegmentTag( int tag ) const
1319 return ( _segTag0 <= tag && tag <= _segTag );
1322 //================================================================================
1324 * \brief Return tag of EDGE by tags of its splits
1326 //================================================================================
1328 int BLSURFPlugin_EnforcedMesh1D::GetTagOfSplitEdge( int splitTag ) const
1330 auto sTag2eTag = _splitTags2EdgeTag.find( splitTag );
1331 return ( sTag2eTag == _splitTags2EdgeTag.end() ) ? splitTag : sTag2eTag->second;