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() ),
151 if ( !hyp || !_mesh || hyp->GetEnforcedMeshes().empty() )
154 _tol = 2 * BRep_Tool::MaxTolerance( _shape, TopAbs_VERTEX );
155 _segTag = 1000 + helper.Count( _shape, TopAbs_EDGE, /*ignoreSame=*/false );
157 _nodeTag0 = 1000 + helper.Count( _shape, TopAbs_VERTEX, /*ignoreSame=*/false );
159 for ( const BLSURFPlugin_Hypothesis::EnforcedMesh& enfMesh : hyp->GetEnforcedMeshes() )
161 copyEnforcedMesh( enfMesh, hyp, _shape );
164 for ( const SMDS_MeshElement* segment : _segmentsOnSeveralShapes )
166 splitSegmentOnSeveralShapes( segment );
169 for ( const TopoDS_Shape & face : _faces )
171 splitSelfIntersectingSegments( face );
174 // for ( TEdge2Nodes::Iterator e2nn( _nodesOnEdge ); e2nn.More(); e2nn.Next() )
176 // splitEdgeByNodes( e2nn.Key(), e2nn.Value() );
180 //================================================================================
182 * \brief Convert enforced segments to quadratic
184 //================================================================================
186 BLSURFPlugin_EnforcedMesh1D::~BLSURFPlugin_EnforcedMesh1D()
191 _helper.SetIsQuadratic( true );
193 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
194 std::vector<const SMDS_MeshNode *> nodes(2);
195 std::vector<const SMDS_MeshElement *> trias;
197 SMDS_MeshGroup* group = nullptr;
198 if ( !_name2Group.empty() )
199 group = _name2Group.begin()->second;
201 for ( const TopoDS_Shape& f : _faces )
202 if ( HasSegmentsOnFace( TopoDS::Face( f )))
203 while ( _segIterator->more() )
205 const SMDS_MeshElement* segment = _segIterator->next();
206 if ( segment->GetType() != SMDSAbs_Edge )
209 nodes[0] = segment->GetNode(0);
210 nodes[1] = segment->GetNode(1);
211 if ( meshDS->GetElementsByNodes( nodes, trias, SMDSAbs_Face ))
213 if ( !group || !group->Contains( segment ))
215 SMDS_MeshGroup* otherGroup = nullptr;
216 if ( _name2Group.size() > bool( group ))
217 for ( auto & n2g : _name2Group )
218 if ( n2g.second != group && n2g.second->Contains( segment ))
220 otherGroup = n2g.second;
226 _helper.AddTLinks( static_cast< const SMDS_MeshFace* >( trias[0] ));
228 meshDS->RemoveFreeElement( segment, meshDS->MeshElements( f ), /*fromGroup=*/false );
230 SMDS_MeshElement* quadSegment = _helper.AddEdge( nodes[0], nodes[1],
231 /*id=*/0, /*force3d=*/false );
232 if ( group && segment != quadSegment )
234 group->Remove( segment );
235 group->Add( quadSegment );
242 //================================================================================
244 * \brief Add a vertex on EDGE
246 //================================================================================
248 void BLSURFPlugin_EnforcedMesh1D::AddVertexOnEdge( const double* theXYZ )
250 // setup predicates to find the supporting EDGE
251 setupPredicates( _shape );
253 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
254 const SMDS_MeshNode* nodeOnE = meshDS->AddNode( theXYZ[0], theXYZ[1], theXYZ[2] );
256 // check if enfNode is on VERTEX
257 bool toRemove = true;
258 TopoDS_Vertex vertex;
260 if ( _onVertexPredicate->IsSatisfy( nodeOnE, &vertex ))
262 toRemove = SMESH_Algo::VertexNode( vertex, meshDS );
264 meshDS->SetNodeOnVertex( nodeOnE, vertex );
266 // find the EDGE supporting theXYZ
267 else if ( _onEdgePredicate->IsSatisfy( nodeOnE, &edge ))
269 gp_Pnt pnt( theXYZ[0], theXYZ[1], theXYZ[2] );
270 toRemove = findNodeOnEdge( pnt, edge );
272 addNodeOnEdge( nodeOnE, edge );
276 meshDS->RemoveFreeNode( nodeOnE, /*submesh=*/nullptr, /*fromGroup=*/false );
279 //================================================================================
281 * \brief Return EDGEs resulted from division of FACE boundary by enforced segments
282 * and enforced vertices
283 * \param [in] edge - possibly divided EDGE
284 * \param [out] splits - split EDGEs
285 * \return bool - true if the EDGE is split
287 //================================================================================
289 bool BLSURFPlugin_EnforcedMesh1D::GetSplitsOfEdge( const TopoDS_Edge& edge,
290 std::vector< TopoDS_Edge > & splits,
291 TopTools_IndexedMapOfShape & edgeTags )
293 if ( _nodesOnEdge.IsBound( edge )) // divide the EDGE
295 splitEdgeByNodes( edge, _nodesOnEdge( edge ));
297 _nodesOnEdge.UnBind( edge );
302 std::vector< TopoDS_Edge > * splitsInMap = _edgeSplitsOfEdge.ChangeSeek( edge );
306 splits.insert( splits.end(), splitsInMap->begin(), splitsInMap->end() );
308 int eTag = edgeTags.Add( edge );
310 size_t index = splits.size() - 1;
311 for ( size_t i = 0; i < splitsInMap->size(); ++i, --index )
313 int splitTag = edgeTags.Add( splits[ index ]);
314 _splitTags2EdgeTag[ splitTag ] = eTag;
316 if ( edge.Orientation() == TopAbs_REVERSED )
317 splits[ index ].Reverse();
324 //================================================================================
326 * \brief Add 1D elements to the mesh
328 //================================================================================
330 void BLSURFPlugin_EnforcedMesh1D::
332 copyEnforcedMesh( const BLSURFPlugin_Hypothesis::EnforcedMesh& theEnfMesh,
333 const BLSURFPlugin_Hypothesis* theHyp,
334 const TopoDS_Shape& theShape)
337 SMDS_ElemIteratorPtr segIt = theHyp->GetEnforcedSegments( theEnfMesh, mesh1D );
338 if ( !segIt->more() )
341 // setup predicates to detect nodes on FACE boundary
342 setupPredicates( theShape );
344 SMDS_MeshGroup* group = getGroup( theEnfMesh._groupName );
345 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
347 // get ordered nodes and segments of theEnfMesh
348 SMESH_MeshAlgos::TElemGroupVector edgeBranches;
349 SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
350 SMESH_MeshAlgos::Get1DBranches( segIt, edgeBranches, nodeBranches );
353 // Copy nodes and segments from an enforced mesh to my mesh
355 TopoDS_Shape vertex, edge;
357 // first treat ends of branches that can be shared by branches
358 for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
360 std::vector< const SMDS_MeshNode* > & braNodes = nodeBranches[ iB ];
361 std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
363 for ( int isLast = 0; isLast < 2; ++isLast )
365 const SMDS_MeshNode* enfNode = isLast ? braNodes.back() : braNodes[0];
366 if ( meshDS->Contains( enfNode ))
367 continue; // already in my mesh
369 const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
371 orientSegments( braSegs, braNodes, isLast ? 0 : braNodes.size() - 1, mesh1D );
373 // replace enfNode at branch ends by newNode
374 SMESH_NodeXYZ enfPnt( newNode ? newNode : enfNode );
375 for ( std::vector< const SMDS_MeshNode* > & braNodes : nodeBranches )
377 for ( int isLast = 0; isLast < 2; ++isLast )
379 const SMDS_MeshNode* & endNode = isLast ? braNodes.back() : braNodes[0];
380 if ( endNode == enfNode || enfPnt.SquareDistance( endNode ) < _tol*_tol )
385 } // loop on branch ends
386 } // loop on nodeBranches
388 // copy nodes and segments
390 for ( size_t iB = 0; iB < nodeBranches.size(); ++iB )
392 std::vector< const SMDS_MeshNode* > & braNodes = nodeBranches[ iB ];
393 std::vector< const SMDS_MeshElement* > & braSegs = edgeBranches[ iB ];
395 // copy nodes of the branch
396 for ( size_t i = 0; i < braNodes.size(); ++i )
398 const SMDS_MeshNode* & enfNode = braNodes[ i ];
399 const SMDS_MeshNode* newNode = copyEnforcedNode( enfNode );
401 if ( !newNode ) // orient segments to be able to get enforced not projected node
402 orientSegments( braSegs, braNodes, i, mesh1D );
407 // copy segments of the branch
408 for ( size_t i = 0; i < braSegs.size(); ++i )
410 //braSegs[ i ] = nullptr;
412 const SMDS_MeshNode* node0 = braNodes[ i ];
413 const SMDS_MeshNode* node1 = braNodes[ i + 1 ];
414 if ( !node0 && !node1 )
417 TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
418 TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
419 if ( shape0.IsNull() && shape1.IsNull() )
422 if ( !node0 && shape1.ShapeType() != TopAbs_FACE )
424 if ( !node1 && shape0.ShapeType() != TopAbs_FACE )
427 if ( !node0 || !node1 ) // create missing node at location of enforced node projected nowhere
429 SMESH_NodeXYZ pn = braSegs[i]->GetNode( !node1 );
430 ( node0 ? node1 : node0 ) = _helper.AddNode( pn->X(), pn->Y(), pn->Z() );
433 SMDS_MeshEdge* newSeg = _helper.AddEdge( node0, node1 );
435 group->Add( newSeg );
436 braSegs[ i ] = newSeg;
438 // check if the both nodes are on the same FACE
439 TopoDS_Shape face = shape0;
440 if ( !shape0.IsSame( shape1 ) && !shape0.IsNull() && !shape1.IsNull() )
442 if ( shape0.ShapeType() == TopAbs_FACE &&
443 _helper.IsSubShape( shape1, shape0 ))
447 else if ( shape1.ShapeType() == TopAbs_FACE &&
448 _helper.IsSubShape( shape0, shape1 ))
452 else // try to find a FACE by projecting a segment middle point
455 gp_Pnt middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
456 //BLSURFPlugin_BLSURF::projectionPoint projPnt =
457 BLSURFPlugin_BLSURF::getProjectionPoint( TopoDS::Face( face ), middlePnt );
459 if ( !face.IsNull() &&
460 ( !_helper.IsSubShape( shape0, face ) ||
461 !_helper.IsSubShape( shape1, face ) ))
465 if ( !face.IsNull() && face.ShapeType() == TopAbs_FACE )
467 meshDS->SetMeshElementOnShape( newSeg, face );
471 if ( face.IsNull() || shape0.IsNull() || shape1.IsNull() )
473 _segmentsOnSeveralShapes.push_back( newSeg );
476 } // loop on branch segments
478 } // loop on branches
483 //================================================================================
485 * \brief Create a copy of a node of enforced mesh in my mesh
486 * \param [in] enfNode - enforced node
487 * \return const SMDS_MeshNode* - a node in my mesh
489 //================================================================================
491 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::copyEnforcedNode( const SMDS_MeshNode* enfNode )
493 SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
495 if ( !enfNode || meshDS->Contains( enfNode ))
496 return enfNode; // already in my mesh
498 SMESH_NodeXYZ enfPnt = enfNode;
500 const SMDS_MeshNode* newNode = nullptr;
502 // check if enfNode is on VERTEX
503 TopoDS_Vertex vertex;
504 if ( _onVertexPredicate->IsSatisfy( enfNode, &vertex ))
506 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
507 newNode = SMESH_Algo::VertexNode( vertex, meshDS );
510 // check if enfNode is on EDGE
511 bool setNodeOnEdge = false;
515 setNodeOnEdge = _onEdgePredicate->IsSatisfy( enfNode, &edge );
518 newNode = findNodeOnEdge( enfPnt, edge );
519 setNodeOnEdge = !newNode;
523 // create a new node and set it on FACE
527 BLSURFPlugin_BLSURF::projectionPoint projPnt =
528 BLSURFPlugin_BLSURF::getProjectionPoint( face, enfPnt, /*allowStateON=*/true );
529 if ( face.IsNull() ) return newNode;
531 if ( projPnt.state == TopAbs_ON )
533 SMDS_MeshNode* projNode = const_cast< SMDS_MeshNode* >( enfNode );
534 projNode->setXYZ( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
537 if ( _onVertexPredicate->IsSatisfy( projNode, &vertex ))
539 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
540 newNode = SMESH_Algo::VertexNode( vertex, meshDS );
542 else if (( setNodeOnEdge = _onEdgePredicate->IsSatisfy( projNode, &edge )))
544 newNode = findNodeOnEdge( projPnt.xyz, edge );
545 setNodeOnEdge = !newNode;
547 projNode->setXYZ( enfPnt.X(), enfPnt.Y(), enfPnt.Z() );
551 newNode = meshDS->AddNode( projPnt.xyz.X(), projPnt.xyz.Y(), projPnt.xyz.Z() );
553 if ( vertex.IsNull() && edge.IsNull() )
554 meshDS->SetNodeOnFace( newNode, face, projPnt.uv.X(), projPnt.uv.Y() );
557 // set the new node on EDGE
558 if ( newNode && setNodeOnEdge )
560 addNodeOnEdge( newNode, edge );
566 //================================================================================
568 * \brief Split a segment whose nodes are on different FACEs into smaller parts
569 * lying each on one FACE
571 //================================================================================
573 void BLSURFPlugin_EnforcedMesh1D::splitSegmentOnSeveralShapes( const SMDS_MeshElement* segment )
575 const SMDS_MeshNode* node0 = segment->GetNode(0);
576 const SMDS_MeshNode* node1 = segment->GetNode(1);
577 SMESHDS_Mesh * meshDS = _helper.GetMeshDS();
579 TopoDS_Shape shape0 = _helper.GetSubShapeByNode( node0, meshDS );
580 TopoDS_Shape shape1 = _helper.GetSubShapeByNode( node1, meshDS );
582 if ( shape0.IsNull() )
584 std::swap( node0, node1 );
589 gp_XYZ xyz0 = SMESH_NodeXYZ( node0 );
590 gp_XYZ xyz1 = SMESH_NodeXYZ( node1 );
591 gp_XYZ segDir = ( xyz1 - xyz0 ).Normalized();
593 //std::map< double, const SMDS_MeshNode* > mediumNodes; // nodes splitting the segment
595 while ( !shape0.IsSame( shape1 )) // move along the segment till shape1
597 if ( shape0.ShapeType() == TopAbs_FACE ) // make a node on an EDGE of the FACE
598 { // ----------------------------------
601 gp_Pnt edgeIntPnt = getEdgeIntersection( shape0, xyz0, xyz1, edge, paramOnE );
605 // check if edgeIntPnt in on VERTEX
606 TopoDS_Vertex vertex;
607 for ( int iV = 0; iV < 2 && vertex.IsNull(); ++iV )
609 TopoDS_Vertex v = _helper.IthVertex( iV, edge );
610 if ( edgeIntPnt.SquareDistance( BRep_Tool::Pnt( v )) < _tol *_tol )
614 // make a node on the EDGE
615 const SMDS_MeshNode* nodeOnEdge = nullptr;
616 if ( vertex.IsNull() )
618 nodeOnEdge = findNodeOnEdge( edgeIntPnt, edge );
621 if ( shape1.IsNull() && node1 )
625 meshDS->MoveNode( nodeOnEdge, edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
629 nodeOnEdge = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
631 addNodeOnEdge( nodeOnEdge, edge, paramOnE );
636 _mesh->GetSubMesh( vertex )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
637 nodeOnEdge = SMESH_Algo::VertexNode( vertex, meshDS );
640 // create a sub-segment
641 SMDS_MeshElement* newSeg = _helper.AddEdge( node0, nodeOnEdge );
642 meshDS->SetMeshElementOnShape( newSeg, shape0 );
644 SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
647 xyz0 = SMESH_NodeXYZ( node0 );
648 if ( vertex.IsNull() )
654 else // shape0 is EDGE or VERTEX; look for the next FACE
655 { // ------------------------------------------------
657 if ( !shape1.IsNull() &&
658 shape1.ShapeType() == TopAbs_FACE &&
659 _helper.IsSubShape( shape0, shape1 )) // shape0 belongs to FACE shape1
661 SMDS_MeshElement* newSeg = _helper.AddEdge( node0, node1 );
662 SMESH_MeshEditor::AddToSameGroups( newSeg, segment, meshDS );
663 meshDS->SetMeshElementOnShape( newSeg, shape1 );
668 double shift = 10 * _tol;
669 for ( int nbAttemp = 0; face.IsNull() && nbAttemp < 10; ++nbAttemp )
671 xyz0 += segDir * shift;
673 BLSURFPlugin_BLSURF::getProjectionPoint( face, xyz0 );
675 if ( !face.IsNull() )
677 if ( _helper.IsSubShape( shape0, face ))
689 } //while ( !shape0.IsSame( shape1 ))
691 meshDS->RemoveFreeElement( segment, /*submesh=*/nullptr );
694 //================================================================================
696 * \brief Find intersection of FACE EDGEs and a segment
697 * \param [in] theFace - the FACE
698 * \param [in] thePnt0 - first end of the segment
699 * \param [in] thePnt1 - last end of the segment
700 * \param [out] theFounfEdge - return the intersected EDGE
701 * \param [out] theParamOnEdge - return parameter of intersection point on EDGE
702 * \return gp_XYZ - point on an EDGE closest to the segment
704 //================================================================================
706 gp_Pnt BLSURFPlugin_EnforcedMesh1D::getEdgeIntersection( const TopoDS_Shape& theFaceOrEdge,
707 const gp_XYZ& thePnt0,
708 const gp_XYZ& thePnt1,
709 TopoDS_Edge & theFounfEdge,
710 double & theParamOnEdge)
712 const double segLen = ( thePnt1 - thePnt0 ).Modulus();
713 const double maxSegDist2 = segLen * segLen * 0.5 * 0.5;
715 Handle(Geom_Line) segLine = new Geom_Line( thePnt0, thePnt1 - thePnt0 );
716 GeomAdaptor_Curve segLineAdpt( segLine, 0, segLen );
718 TopTools_MapOfShape edges;
719 double minParamOnSeg = segLen;
721 for ( TopExp_Explorer edgeExp( theFaceOrEdge, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
723 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
724 if ( !edges.Add( edge ))
727 Extrema_ExtCC extrema( segLineAdpt, BRepAdaptor_Curve( edge ), _tol, _tol );
729 if ( extrema.IsDone() && !extrema.IsParallel() )
730 for ( int i = 1, nb = extrema.NbExt(); i <= nb; ++i )
731 if ( extrema.SquareDistance( i ) < maxSegDist2 )
733 Extrema_POnCurv pOnSeg, pOnEdge;
734 extrema.Points( i, pOnSeg, pOnEdge );
735 double paramOnSeg = pOnSeg.Parameter();
736 if ( 0 < paramOnSeg && paramOnSeg < minParamOnSeg )
738 minParamOnSeg = paramOnSeg;
739 foundPnt = pOnEdge.Value();
741 theParamOnEdge = pOnEdge.Parameter();
748 //================================================================================
750 * \brief Split self-intersecting segments on a given FACE
752 //================================================================================
754 void BLSURFPlugin_EnforcedMesh1D::splitSelfIntersectingSegments( const TopoDS_Shape & theFace )
756 const TopoDS_Face& face = TopoDS::Face( theFace );
758 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
759 SMESHDS_SubMesh* sm = meshDS->MeshElements( face );
760 if ( !sm || sm->NbElements() <= 1 )
763 // get ordered nodes and segments on the face
764 SMESH_MeshAlgos::TElemGroupVector edgeBranches;
765 SMESH_MeshAlgos::TNodeGroupVector nodeBranches;
766 SMESH_MeshAlgos::Get1DBranches( sm->GetElements(), edgeBranches, nodeBranches );
768 // create element searcher
769 SMESH_ElementSearcher* elemSearcher;
770 SMESHUtils::Deleter< SMESH_ElementSearcher > elSearchdeleter;
771 std::vector< const SMDS_MeshElement* > foundElems;
773 std::vector< SMDS_ElemIteratorPtr > elemItVec;
774 for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
776 SMDS_ElemIteratorPtr segIt =
777 boost::make_shared< SMDS_ElementVectorIterator >( braSegs.begin(), braSegs.end() );
778 elemItVec.push_back( segIt );
780 typedef SMDS_IteratorOnIterators< const SMDS_MeshElement*,
781 std::vector< SMDS_ElemIteratorPtr > > TVecIterator;
782 SMDS_ElemIteratorPtr segIt = boost::make_shared< TVecIterator >( elemItVec );
784 elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, segIt, _tol );
785 elSearchdeleter._obj = elemSearcher;
787 // force usage of iterators before they die
788 elemSearcher->FindElementsByPoint( gp_Pnt( 0,0,1e+20), SMDSAbs_Edge, foundElems );
792 // Find intersecting segments
794 std::map< const SMDS_MeshElement* , std::vector< const SMDS_MeshNode* > > segInternalNodes;
795 SMESH_MeshEditor::TListOfListOfNodes nodeGroupsToMerge;
797 for ( std::vector< const SMDS_MeshElement* > & braSegs : edgeBranches )
799 braSegs.push_back( braSegs.back() );
800 const SMDS_MeshElement* prevSeg = nullptr;
802 for ( size_t i = 0, nb = braSegs.size() - 1; i < nb; ++i )
804 const SMDS_MeshElement* seg = braSegs[ i ];
805 const SMDS_MeshElement* nextSeg = braSegs[ i + 1 ];
806 const SMDS_MeshNode* node0 = seg->GetNode(0);
807 const SMDS_MeshNode* node1 = seg->GetNode(1);
809 gp_XYZ xyz0 = SMESH_NodeXYZ( node0 );
810 gp_XYZ xyz1 = SMESH_NodeXYZ( node1 );
811 gp_XYZ middlePnt = 0.5 * ( SMESH_NodeXYZ( node0 ) + SMESH_NodeXYZ( node1 ));
812 double segLen = ( xyz0 - xyz1 ).Modulus();
815 elemSearcher->GetElementsInSphere( middlePnt, 0.5 * segLen + _tol, SMDSAbs_Edge, foundElems );
817 for ( const SMDS_MeshElement* closeSeg : foundElems )
819 if ( closeSeg == prevSeg ||
821 closeSeg == nextSeg )
826 if ( intersectSegments( seg, closeSeg, face, intPnt, uv ))
829 const SMDS_MeshNode* intNode0 = findNode( intPnt, _tol, seg, segInternalNodes, i0 );
830 const SMDS_MeshNode* intNode1 = findNode( intPnt, _tol, closeSeg, segInternalNodes, i1 );
831 if ( !intNode0 && intNode1 )
833 segInternalNodes[ seg ].push_back( intNode1 );
835 else if ( !intNode1 && intNode0 )
837 segInternalNodes[ closeSeg ].push_back( intNode0 );
839 else if ( intNode1 && intNode0 )
841 if ( intNode1 == intNode0 )
843 if ( i0 < 0 && i1 < 0 )
845 nodeGroupsToMerge.push_back // merge end nodes
846 ( std::list< const SMDS_MeshNode* >({ intNode0, intNode1 }));
850 segInternalNodes[ closeSeg ][ i1 ] = intNode0;
854 segInternalNodes[ seg ][ i0 ] = intNode1;
856 else // two internal nodes coincide
858 segInternalNodes[ seg ][ i0 ] = intNode1;
859 nodeGroupsToMerge.push_back
860 ( std::list< const SMDS_MeshNode* >({ intNode1, intNode0 }));
863 else // ( !intNode1 && !intNode0 )
865 intNode0 = _helper.AddNode( intPnt.X(), intPnt.Y(), intPnt.Z() );
866 meshDS->SetNodeOnFace( intNode0, face, uv.X(), uv.Y() );
867 segInternalNodes[ seg ].push_back( intNode0 );
868 segInternalNodes[ closeSeg ].push_back( intNode0 );
875 } // loop on segments of a branch
876 } // loop on branches
879 findIntersectionWithSeamEdge( face, segInternalNodes ); // on periodic FACE
884 for ( auto& seg2nodes : segInternalNodes )
886 const SMDS_MeshElement* seg = seg2nodes.first;
887 std::vector< const SMDS_MeshNode* > & nodes = seg2nodes.second;
888 if ( nodes.empty() ) continue;
890 const SMDS_MeshNode* n0 = seg->GetNode( 0 );
891 const SMDS_MeshNode* n1 = seg->GetNode( 1 );
892 nodes.push_back( n1 );
894 // sort nodes on the segment
895 gp_Pnt p0 = SMESH_NodeXYZ( n0 );
896 std::map< double, const SMDS_MeshNode* > sortedNodes;
897 for ( SMESH_NodeXYZ pn : nodes )
898 sortedNodes.insert({ p0.SquareDistance( pn ), pn.Node() });
901 for ( auto & d2n : sortedNodes )
904 SMDS_MeshElement* newSeg = _helper.AddEdge( n0, n1 );
906 meshDS->SetMeshElementOnShape( newSeg, face );
907 SMESH_MeshEditor::AddToSameGroups( newSeg, seg, meshDS );
909 meshDS->RemoveFreeElement( seg, /*submesh=*/nullptr );
914 SMESH_MeshEditor( _mesh ).MergeNodes( nodeGroupsToMerge );
917 //================================================================================
919 * \brief Find intersections of segments with a seam EDGE on a periodic FACE
921 //================================================================================
923 void BLSURFPlugin_EnforcedMesh1D::findIntersectionWithSeamEdge( const TopoDS_Face & face,
924 TNodesOnSeg & segInternalNodes )
926 _helper.SetSubShape( face );
927 if ( !_helper.HasSeam() )
930 SMESHDS_Mesh* meshDS = _helper.GetMeshDS();
931 SMESHDS_SubMesh* sm = meshDS->MeshElements( face );
932 if ( !sm || sm->NbElements() == 0 )
935 TopTools_MapOfShape treatedEdges;
936 for ( TopExp_Explorer edgeExp( face, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
938 const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
939 if ( !_helper.IsSeamShape( edge ) || !treatedEdges.Add( edge ))
942 for ( SMDS_ElemIteratorPtr setIt = sm->GetElements(); setIt->more(); )
944 const SMDS_MeshElement* segment = setIt->next();
945 const SMESH_NodeXYZ pn0 = segment->GetNode( 0 );
946 const SMESH_NodeXYZ pn1 = segment->GetNode( 1 );
948 gp_XY uv0 = _helper.GetNodeUV( face, pn0.Node() );
949 gp_XY uv1 = _helper.GetNodeUV( face, pn1.Node() );
951 for ( int iCoo = 1; iCoo <= 2; ++iCoo )
952 if ( iCoo & _helper.GetPeriodicIndex() )
954 double distParam = Abs( uv0.Coord( iCoo ) - uv1.Coord( iCoo ));
955 if ( distParam > 0.5 * _helper.GetPeriod( iCoo ))
957 double paramOnE; TopoDS_Edge edge2;
958 gp_Pnt edgeIntPnt = getEdgeIntersection( edge, pn0, pn1, edge2, paramOnE );
959 if ( edge2.IsNull() )
962 const SMDS_MeshNode* nodeOnSeam = findNodeOnEdge( edgeIntPnt, edge );
963 bool isOnEdge = nodeOnSeam;
965 int indexOnSegment = 3;
967 nodeOnSeam = findNode( edgeIntPnt, _tol, segment, segInternalNodes, indexOnSegment );
971 nodeOnSeam = _helper.AddNode( edgeIntPnt.X(), edgeIntPnt.Y(), edgeIntPnt.Z() );
976 addNodeOnEdge( nodeOnSeam, edge, paramOnE );
978 if ( indexOnSegment == 3 )
979 segInternalNodes[ segment ].push_back( nodeOnSeam );
986 //================================================================================
988 * \brief Intersect two segments
989 * \param [in] seg1 - segment 1
990 * \param [in] seg2 - segment 2
991 * \param [in] face - the FACE on which segments lie
992 * \param [out] intPnt - intersection point
993 * \param [out] intUV - UV of intersection point on the FACE
994 * \return bool - true if intersection found
996 //================================================================================
998 bool BLSURFPlugin_EnforcedMesh1D::intersectSegments( const SMDS_MeshElement* seg1,
999 const SMDS_MeshElement* seg2,
1000 const TopoDS_Face& face,
1002 gp_Pnt2d & intUV ) const
1004 SMESH_NodeXYZ n10 = seg1->GetNode(0);
1005 SMESH_NodeXYZ n11 = seg1->GetNode(1);
1006 SMESH_NodeXYZ n20 = seg2->GetNode(0);
1007 SMESH_NodeXYZ n21 = seg2->GetNode(1);
1008 if ( n10 == n20 || n10 == n21 || n11 == n20 || n10 == n21 )
1011 gp_Lin lin1( n10, n11 - n10 );
1012 gp_Lin lin2( n20, n21 - n20 );
1014 Extrema_ExtElC extrema( lin1, lin2, Precision::Angular() );
1016 if ( !extrema.IsDone() || extrema.IsParallel() )
1019 Extrema_POnCurv poc1, poc2;
1020 extrema.Points( 1, poc1, poc2 );
1022 double len1 = lin1.Direction().XYZ() * ( n11 - n10 );
1023 double len2 = lin2.Direction().XYZ() * ( n21 - n20 );
1024 double u1 = poc1.Parameter();
1025 double u2 = poc2.Parameter();
1027 if ( u1 < -_tol || u1 > len1 + _tol )
1029 if ( u2 < -_tol || u2 > len2 + _tol )
1032 intPnt = 0.5 * ( poc1.Value().XYZ() + poc2.Value().XYZ() );
1034 // compute approximate UV
1039 gp_XY uv10 = _helper.GetNodeUV( face, n10.Node(), n11.Node() );
1040 gp_XY uv11 = _helper.GetNodeUV( face, n11.Node(), n10.Node() );
1041 gp_XY uv1 = uv10 * ( 1 - u1 ) + uv11 * u1;
1043 gp_XY uv20 = _helper.GetNodeUV( face, n20.Node(), n21.Node() );
1044 gp_XY uv21 = _helper.GetNodeUV( face, n21.Node(), n20.Node() );
1045 gp_XY uv2 = uv20 * ( 1 - u2 ) + uv21 * u2;
1047 intUV = 0.5 * ( uv1 + uv2 );
1049 // compute precise UV and XYZ by projecting intPnt to the FACE
1051 Handle(ShapeAnalysis_Surface) surface = _helper.GetSurface( face );
1052 intUV = surface->NextValueOfUV( intUV, intPnt, _tol );
1053 if ( surface->Gap() > Min( len1, len2 ))
1054 intUV = surface->ValueOfUV( intPnt, _tol );
1056 intPnt = surface->Value( intUV );
1061 //================================================================================
1063 * \brief Setup predicates to detect nodes on FACE boundary
1064 * \param [in] shape - shape containing FACEs to mesh
1066 //================================================================================
1068 void BLSURFPlugin_EnforcedMesh1D::setupPredicates( const TopoDS_Shape& theShape )
1070 if ( _onVertexPredicate )
1073 _onVertexPredicate.reset( new TPredicate() );
1074 _onVertexPredicate->SetTolerance( _tol );
1075 _onVertexPredicate->SetMesh( _helper.GetMeshDS() );
1077 TopTools_IndexedMapOfShape vertices;
1078 TopExp::MapShapes( theShape, TopAbs_VERTEX, vertices );
1079 TopoDS_Compound vCompound;
1080 BRep_Builder builder;
1081 builder.MakeCompound( vCompound );
1082 for ( const TopoDS_Shape& v : vertices )
1083 builder.Add( vCompound, v );
1085 _onVertexPredicate->SetShape( vCompound, SMDSAbs_Node );
1089 _onEdgePredicate.reset( new TPredicate() );
1090 _onEdgePredicate->SetTolerance( _tol );
1091 _onEdgePredicate->SetMesh( _helper.GetMeshDS() );
1093 TopTools_IndexedMapOfShape edges;
1094 TopExp::MapShapes( theShape, TopAbs_EDGE, edges );
1095 TopoDS_Compound eCompound;
1096 BRep_Builder builder;
1097 builder.MakeCompound( eCompound );
1098 for ( const TopoDS_Shape& e : edges )
1099 builder.Add( eCompound, e );
1101 _onEdgePredicate->SetShape( eCompound, SMDSAbs_Node );
1106 //================================================================================
1108 * \brief Find or create a group of edges with given name
1110 //================================================================================
1112 SMDS_MeshGroup* BLSURFPlugin_EnforcedMesh1D::getGroup( const std::string& groupName )
1114 SMDS_MeshGroup* group = nullptr;
1115 if ( !groupName.empty() )
1117 if ( _name2Group.count( groupName ))
1118 return _name2Group[ groupName ];
1120 // find existing group
1121 for ( SMESH_Mesh::GroupIteratorPtr grIt = _mesh->GetGroups(); grIt->more(); )
1123 SMESH_Group* grp = grIt->next();
1124 SMESHDS_Group* grpDS = dynamic_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1126 grpDS->GetType() == SMDSAbs_Edge &&
1127 groupName == grp->GetName() )
1129 _name2Group[ groupName ] = & grpDS->SMDSGroup();
1130 return & grpDS->SMDSGroup();
1134 // create a new group
1135 SMESH_Group* grp = _mesh->AddGroup( SMDSAbs_Edge, groupName.c_str() );
1136 SMESHDS_Group* grpDS = static_cast< SMESHDS_Group* >( grp->GetGroupDS() );
1138 group = & grpDS->SMDSGroup();
1139 _name2Group[ groupName ] = group;
1144 //================================================================================
1146 * \brief Look for a node dividing a given EDGE
1148 //================================================================================
1150 const SMDS_MeshNode* BLSURFPlugin_EnforcedMesh1D::findNodeOnEdge( const gp_Pnt& p,
1151 const TopoDS_Edge& edge )
1153 // look for an equal node on the EDGE
1154 if ( std::vector< const SMDS_MeshNode* >* nodesOnE = _nodesOnEdge.ChangeSeek( edge ))
1155 for ( const SMDS_MeshNode* n : *nodesOnE )
1156 if ( p.SquareDistance( SMESH_NodeXYZ( n )) < _tol * _tol )
1164 //================================================================================
1166 * \brief Add a node to an EDGE
1168 //================================================================================
1170 void BLSURFPlugin_EnforcedMesh1D::addNodeOnEdge( const SMDS_MeshNode* node,
1171 const TopoDS_Edge& edge,
1174 _mesh->GetMeshDS()->SetNodeOnEdge( node, edge, u );
1176 std::vector< const SMDS_MeshNode* > * nodesOnE = _nodesOnEdge.ChangeSeek( edge );
1178 nodesOnE = _nodesOnEdge.Bound( edge, std::vector< const SMDS_MeshNode* >() );
1180 nodesOnE->push_back( node );
1183 //================================================================================
1185 * \brief Create EDGEs by dividing a given EDGE by nodes on it
1186 * \param [in] edge - the EDGE to divide
1187 * \param [in] nodes - the nodes to divide by
1189 //================================================================================
1191 void BLSURFPlugin_EnforcedMesh1D::
1193 splitEdgeByNodes( TopoDS_Edge edge, const std::vector< const SMDS_MeshNode* >& nodes )
1195 if ( nodes.empty() )
1198 edge.Orientation( TopAbs_FORWARD );
1200 TopoDS_Vertex v0 = _helper.IthVertex( 0, edge );
1201 TopoDS_Vertex v1 = _helper.IthVertex( 1, edge );
1203 // create VERTEXes and sort them along the EDGE
1205 std::map< double, TopoDS_Vertex > sortedVertices;
1207 BRepAdaptor_Curve curve( edge );
1208 for ( SMESH_NodeXYZ pn : nodes )
1212 ShapeAnalysis_Curve().Project( curve, pn, _tol, projPnt, u, false );
1213 projPnt = curve.Value( u );
1215 TopoDS_Vertex v = BRepBuilderAPI_MakeVertex( projPnt );
1217 sortedVertices.insert({ u, v });
1219 _nOnE2Vertex[ pn.Node() ] = v;
1221 sortedVertices.insert({ BRep_Tool::Parameter( v1, edge ), v1 });
1226 BRep_Builder builder;
1227 std::vector< TopoDS_Edge >& newEdges = *_edgeSplitsOfEdge.Bound( edge, TEdge2Edges::value_type());
1229 double u0 = BRep_Tool::Parameter( v0, edge );
1230 for ( auto& u2v : sortedVertices )
1232 double u1 = u2v.first;
1235 TopoDS_Shape newShape = edge.EmptyCopied();
1236 TopoDS_Edge newEdge = TopoDS::Edge( newShape );
1237 builder.Add( newEdge, v0 );
1238 builder.Add( newEdge, v1 );
1239 builder.Range( newEdge, u0, u1 );
1240 newEdges.push_back( newEdge );
1249 //================================================================================
1251 * \brief Return true if there are enforced segments on a FACE. Start iteration on segments
1253 //================================================================================
1255 bool BLSURFPlugin_EnforcedMesh1D::HasSegmentsOnFace( const TopoDS_Face& face )
1257 _segIterator.reset();
1259 if ( _faces.Contains( face ))
1261 _currentFace = face;
1262 _segIterator = _helper.GetMeshDS()->MeshElements( face )->GetElements();
1264 _helper.SetSubShape( face );
1266 return _segIterator->more();
1272 //================================================================================
1274 * \brief Return next segment on the FACE
1276 //================================================================================
1278 bool BLSURFPlugin_EnforcedMesh1D::NextSegment( Segmemnt & seg,
1279 TopTools_IndexedMapOfShape & vertexTags )
1281 if ( _segIterator && _segIterator->more() )
1283 const SMDS_MeshElement* segment = _segIterator->next();
1285 while ( segment->GetType() != SMDSAbs_Edge && _segIterator->more() )
1286 segment = _segIterator->next();
1287 if ( segment->GetType() != SMDSAbs_Edge )
1290 const SMDS_MeshNode * node[2] = { segment->GetNode(0),
1291 segment->GetNode(1) };
1293 seg._tag = _segTag++;
1295 seg._xyz[0] = SMESH_NodeXYZ( node[0]);
1296 seg._xyz[1] = SMESH_NodeXYZ( node[1]);
1298 seg._uv[0] = _helper.GetNodeUV( _currentFace, node[0], node[1] );
1299 seg._uv[1] = _helper.GetNodeUV( _currentFace, node[1], node[0] );
1301 seg._pcurve = new Geom2d_Line( seg._uv[0], ( seg._uv[1] - seg._uv[0] ));
1304 seg._u[1] = ( seg._uv[1] - seg._uv[0] ).Modulus();
1306 for ( int i = 0; i < 2; ++i )
1308 auto n2v = _nOnE2Vertex.find( node[i] ); // find VERTEX by node
1310 if ( n2v != _nOnE2Vertex.end() )
1311 seg._vTag[i] = vertexTags.Add( n2v->second );
1313 seg._vTag[i] = _nodeTag0 + _nodeTags.Add( node[i] );
1316 if ( !_isQuadratic )
1317 _mesh->GetMeshDS()->UnSetMeshElementOnShape( segment, _currentFace );
1324 //================================================================================
1326 * \brief Return enforced node by the tag that was returned by Segmemnt::_vTag[i]
1328 //================================================================================
1330 const SMDS_MeshNode*
1331 BLSURFPlugin_EnforcedMesh1D::GetNodeByTag( int tag,
1332 const TopTools_IndexedMapOfShape & vertexTags )
1334 const SMDS_MeshNode* node = nullptr;
1336 if ( tag <= vertexTags.Size() && !_nOnE2Vertex.empty() )
1338 const TopoDS_Shape& vertex = vertexTags( tag );
1340 for ( auto n2v = _nOnE2Vertex.begin(); n2v != _nOnE2Vertex.end(); ++n2v )
1341 if ( vertex.IsSame( n2v->second ))
1344 _nOnE2Vertex.erase( n2v );
1351 bool isValid = ( 0 < tag && tag <= _nodeTags.Size() );
1353 node = _nodeTags( tag );
1358 //================================================================================
1360 * \brief Return true if a tag corresponds to the tag of enforced segment
1361 * that was returned by Segment::_tag
1363 //================================================================================
1365 bool BLSURFPlugin_EnforcedMesh1D::IsSegmentTag( int tag ) const
1367 return ( _segTag0 <= tag && tag <= _segTag );
1370 //================================================================================
1372 * \brief Return tag of EDGE by tags of its splits
1374 //================================================================================
1376 int BLSURFPlugin_EnforcedMesh1D::GetTagOfSplitEdge( int splitTag ) const
1378 auto sTag2eTag = _splitTags2EdgeTag.find( splitTag );
1379 return ( sTag2eTag == _splitTags2EdgeTag.end() ) ? splitTag : sTag2eTag->second;