1 // Copyright (C) 2018-2019 OPEN CASCADE
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 : SMESH_Slot.cxx
20 // Created : Fri Nov 30 15:58:37 2018
21 // Author : Edward AGAPOV (eap)
23 #include "SMESH_MeshAlgos.hxx"
25 #include "ObjectPool.hxx"
26 #include "SMDS_LinearEdge.hxx"
27 #include "SMDS_Mesh.hxx"
28 #include "SMDS_MeshGroup.hxx"
30 #include <IntAna_IntConicQuad.hxx>
31 #include <IntAna_Quadric.hxx>
32 #include <NCollection_DataMap.hxx>
33 #include <NCollection_Map.hxx>
34 #include <Precision.hxx>
36 #include <gp_Cylinder.hxx>
43 #include <Utils_SALOME_Exception.hxx>
47 typedef SMESH_MeshAlgos::Edge TEdge;
49 //================================================================================
50 //! point of intersection of a face edge with the cylinder
53 SMESH_NodeXYZ myNode; // point and a node
54 int myEdgeIndex; // face edge index
55 bool myIsOutPln[2]; // isOut of two planes
57 double SquareDistance( const IntPoint& p ) const { return ( myNode-p.myNode ).SquareModulus(); }
60 //================================================================================
64 IntPoint myIntPnt1, myIntPnt2;
65 const SMDS_MeshElement* myFace;
67 const IntPoint& operator[]( size_t i ) const { return i ? myIntPnt2 : myIntPnt1; }
69 double SquareDistance( const gp_Pnt& p, gp_XYZ & pClosest ) const
71 gp_Vec edge( myIntPnt1.myNode, myIntPnt2.myNode );
72 gp_Vec n1p ( myIntPnt1.myNode, p );
73 double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
76 pClosest = myIntPnt1.myNode;
77 return n1p.SquareMagnitude();
81 pClosest = myIntPnt2.myNode;
82 return p.SquareDistance( myIntPnt2.myNode );
84 pClosest = myIntPnt1.myNode + u * edge.XYZ(); // projection of the point on the edge
85 return p.SquareDistance( pClosest );
89 //================================================================================
93 typedef std::vector< Cut > TCutList;
95 const SMDS_MeshElement* myEdge;
97 std::vector< const IntPoint* > myFreeEnds; // ends of cut edges
99 Segment( const SMDS_MeshElement* e = 0 ): myEdge(e) { myCuts.reserve( 4 ); }
102 gp_Ax1 Ax1( bool reversed = false ) const
104 SMESH_NodeXYZ n1 = myEdge->GetNode( reversed );
105 SMESH_NodeXYZ n2 = myEdge->GetNode( !reversed );
106 return gp_Ax1( n1, gp_Dir( n2 - n1 ));
110 const SMDS_MeshNode* Node(int i) const
112 return myEdge->GetNode( i % 2 );
115 // store an intersection edge forming the slot border
116 void AddCutEdge( const IntPoint& p1,
118 const SMDS_MeshElement* myFace )
120 myCuts.push_back( Cut({ p1, p2, myFace }));
123 // return number of not shared IntPoint's
124 int NbFreeEnds( double tol )
126 if ( myCuts.empty() )
128 if ( myFreeEnds.empty() )
130 // remove degenerated cuts
131 // for ( size_t iC1 = 0; iC1 < myCuts.size(); ++iC1 )
132 // if ( myCuts[ iC1 ][ 0 ].myNode == myCuts[ iC1 ][ 1 ].myNode )
134 // if ( iC1 < myCuts.size() - 1 )
135 // myCuts[ iC1 ] = myCuts.back();
136 // myCuts.pop_back();
140 std::vector< bool > isSharedPnt( myCuts.size() * 2, false );
141 for ( size_t iC1 = 0; iC1 < myCuts.size() - 1; ++iC1 )
142 for ( size_t iP1 = 0; iP1 < 2; ++iP1 )
144 size_t i1 = iC1 * 2 + iP1;
145 if ( isSharedPnt[ i1 ])
147 for ( size_t iC2 = iC1 + 1; iC2 < myCuts.size(); ++iC2 )
148 for ( size_t iP2 = 0; iP2 < 2; ++iP2 )
150 size_t i2 = iC2 * 2 + iP2;
151 if ( isSharedPnt[ i2 ])
153 if ( myCuts[ iC1 ][ iP1 ].SquareDistance( myCuts[ iC2 ][ iP2 ]) < tol * tol )
156 isSharedPnt[ i1 ] = isSharedPnt[ i2 ] = true;
160 myFreeEnds.reserve( isSharedPnt.size() - nbShared );
161 for ( size_t i = 0; i < isSharedPnt.size(); ++i )
162 if ( !isSharedPnt[ i ] )
166 myFreeEnds.push_back( & myCuts[ iC ][ iP ]);
169 return myFreeEnds.size();
172 typedef ObjectPoolIterator<Segment> TSegmentIterator;
175 //================================================================================
176 //! Segments and plane separating domains of segments, at common node
179 std::vector< Segment* > mySegments;
180 gp_Ax1 myPlane; // oriented OK for mySegments[0]
182 void AddSegment( Segment* seg, const SMDS_MeshNode* n )
184 mySegments.reserve(2);
185 mySegments.push_back( seg );
186 if ( mySegments.size() == 1 )
188 myPlane = mySegments[0]->Ax1( mySegments[0]->myEdge->GetNodeIndex( n ));
192 gp_Ax1 axis2 = mySegments[1]->Ax1( mySegments[1]->myEdge->GetNodeIndex( n ));
193 myPlane.SetDirection( myPlane.Direction().XYZ() - axis2.Direction().XYZ() );
196 gp_Ax1 Plane( const Segment* seg )
198 return ( seg == mySegments[0] ) ? myPlane : myPlane.Reversed();
201 typedef NCollection_DataMap< const SMDS_MeshNode*, NodeData, SMESH_Hasher > TSegmentsOfNode;
204 //================================================================================
206 * \brief Intersect a face edge given by its nodes with a cylinder.
208 //================================================================================
210 bool intersectEdge( const gp_Cylinder& cyl,
211 const SMESH_NodeXYZ& n1,
212 const SMESH_NodeXYZ& n2,
214 std::vector< IntPoint >& intPoints )
216 gp_Lin line( gp_Ax1( n1, gp_Dir( n2 - n1 )));
217 IntAna_IntConicQuad intersection( line, IntAna_Quadric( cyl ));
219 if ( !intersection.IsDone() ||
220 intersection.IsParallel() ||
221 intersection.IsInQuadric() ||
222 intersection.NbPoints() == 0 )
225 gp_Vec edge( n1, n2 );
227 size_t oldNbPnts = intPoints.size();
228 for ( int iP = 1; iP <= intersection.NbPoints(); ++iP )
230 const gp_Pnt& p = intersection.Point( iP );
232 gp_Vec n1p ( n1, p );
233 const SMDS_MeshNode* n = 0;
235 double u = ( edge * n1p ) / edge.SquareMagnitude(); // param [0,1] on the edge
237 if ( p.SquareDistance( n1 ) < tol * tol )
242 else if ( u >= 1. ) {
243 if ( p.SquareDistance( n2 ) < tol * tol )
249 if ( p.SquareDistance( n1 ) < tol * tol )
251 else if ( p.SquareDistance( n2 ) < tol * tol )
255 intPoints.push_back( IntPoint() );
257 intPoints.back().myNode.Set( n );
259 intPoints.back().myNode.SetCoord( p.X(),p.Y(),p.Z() );
262 // set points order along an edge
263 if ( intPoints.size() - oldNbPnts == 2 &&
264 intersection.ParamOnConic( 1 ) > intersection.ParamOnConic( 2 ))
266 int i = intPoints.size() - 1;
267 std::swap( intPoints[ i ], intPoints[ i - 1 ]);
270 return intPoints.size() - oldNbPnts > 0;
273 //================================================================================
275 * \brief Return signed distance between a point and a plane
277 //================================================================================
279 double signedDist( const gp_Pnt& p, const gp_Ax1& planeNormal )
281 const gp_Pnt& O = planeNormal.Location();
283 return Op * planeNormal.Direction();
286 //================================================================================
288 * \brief Check if a point is outside a segment domain bound by two planes
290 //================================================================================
292 bool isOut( const gp_Pnt& p, const gp_Ax1* planeNormal, bool* isOutPtr, int nbPln = 2 )
294 isOutPtr[0] = isOutPtr[1] = false;
296 for ( int i = 0; i < nbPln; ++i )
298 isOutPtr[i] = ( signedDist( p, planeNormal[i] ) <= 0. );
300 return ( isOutPtr[0] && isOutPtr[1] );
303 //================================================================================
305 * \brief Check if a segment between two points is outside a segment domain bound by two planes
307 //================================================================================
309 bool isSegmentOut( bool* isOutPtr1, bool* isOutPtr2 )
311 return (( isOutPtr1[0] && isOutPtr2[0] ) ||
312 ( isOutPtr1[1] && isOutPtr2[1] ));
315 //================================================================================
317 * \brief cut off ip1 from edge (ip1 - ip2) by a plane
319 //================================================================================
321 void cutOff( IntPoint & ip1, const IntPoint & ip2, const gp_Ax1& planeNormal, double tol )
323 gp_Lin lin( ip1.myNode, ( ip2.myNode - ip1.myNode ));
324 gp_Pln pln( planeNormal.Location(), planeNormal.Direction() );
326 IntAna_IntConicQuad intersection( lin, pln, Precision::Angular/*Tolerance*/() );
327 if ( intersection.IsDone() &&
328 !intersection.IsParallel() &&
329 !intersection.IsInQuadric() &&
330 intersection.NbPoints() == 1 )
332 if ( intersection.Point( 1 ).SquareDistance( ip1.myNode ) > tol * tol )
334 static_cast< gp_XYZ& >( ip1.myNode ) = intersection.Point( 1 ).XYZ();
335 ip1.myNode._node = 0;
336 ip1.myEdgeIndex = -1;
341 //================================================================================
343 * \brief Assure that face normal is computed in faceNormals vector
345 //================================================================================
347 const gp_XYZ& computeNormal( const SMDS_MeshElement* face,
348 std::vector< gp_XYZ >& faceNormals )
351 if ((int) faceNormals.size() <= face->GetID() )
354 faceNormals.resize( face->GetID() + 1 );
358 toCompute = faceNormals[ face->GetID() ].SquareModulus() == 0.;
361 SMESH_MeshAlgos::FaceNormal( face, faceNormals[ face->GetID() ], /*normalized=*/false );
363 return faceNormals[ face->GetID() ];
366 typedef std::vector< SMDS_MeshGroup* > TGroupVec;
368 //================================================================================
370 * \brief Fill theFaceID2Groups map for a given face
371 * \param [in] theFace - the face
372 * \param [in] theGroupsToUpdate - list of groups to treat
373 * \param [out] theFaceID2Groups - the map to fill in
374 * \param [out] theWorkGroups - a working buffer of groups
376 //================================================================================
378 void findGroups( const SMDS_MeshElement * theFace,
379 TGroupVec & theGroupsToUpdate,
380 NCollection_DataMap< int, TGroupVec > & theFaceID2Groups,
381 TGroupVec & theWorkGroups )
383 theWorkGroups.clear();
384 for ( size_t i = 0; i < theGroupsToUpdate.size(); ++i )
385 if ( theGroupsToUpdate[i]->Contains( theFace ))
386 theWorkGroups.push_back( theGroupsToUpdate[i] );
388 if ( !theWorkGroups.empty() )
389 theFaceID2Groups.Bind( theFace->GetID(), theWorkGroups );
392 //================================================================================
394 * \brief Check distance between a point and an edge defined by a couple of nodes
396 //================================================================================
398 bool isOnEdge( const SMDS_MeshNode* n1,
399 const SMDS_MeshNode* n2,
403 SMDS_LinearEdge edge( n1, n2 );
404 return ( SMESH_MeshAlgos::GetDistance( &edge, p ) < tol );
407 //================================================================================
409 * \return Index of intersection point detected on a triangle cut by planes
410 * \param [in] i - index of a cut triangle side
411 * \param [in] n1 - 1st point of a cut triangle side
412 * \param [in] n2 - 2nd point of a cut triangle side
413 * \param [in] face - a not cut triangle
414 * \param [in] intPoint - the intersection point
415 * \param [in] faceNodes - nodes of not cut triangle
416 * \param [in] tol - tolerance
418 //================================================================================
420 int edgeIndex( const int i,
421 const SMESH_NodeXYZ& n1,
422 const SMESH_NodeXYZ& n2,
423 const SMDS_MeshElement* face,
424 const IntPoint& intPoint,
425 const std::vector< const SMDS_MeshNode* >& faceNodes,
428 if ( n1.Node() && n2.Node() )
429 return face->GetNodeIndex( n1.Node() );
431 // project intPoint to sides of face
432 for ( size_t i = 1; i < faceNodes.size(); ++i )
433 if ( isOnEdge( faceNodes[ i-1 ], faceNodes[ i ], intPoint.myNode, tol ))
439 //================================================================================
441 * \brief Find a neighboring segment and its next node
442 * \param [in] curSegment - a current segment
443 * \param [in,out] curNode - a current node to update
444 * \param [in] segmentsOfNode - map of segments of nodes
445 * \return Segment* - the found segment
447 //================================================================================
449 Segment* nextSegment( const Segment* curSegment,
450 const SMDS_MeshNode* & curNode,
451 const TSegmentsOfNode& segmentsOfNode )
453 Segment* neighborSeg = 0;
454 const NodeData& noData = segmentsOfNode( curNode );
455 for ( size_t iS = 0; iS < noData.mySegments.size() && !neighborSeg; ++iS )
456 if ( noData.mySegments[ iS ] != curSegment )
457 neighborSeg = noData.mySegments[ iS ];
461 int iN = ( neighborSeg->Node(0) == curNode );
462 curNode = neighborSeg->Node( iN );
467 //================================================================================
469 * \brief Tries to find a segment to which a given point is too close
470 * \param [in] p - the point
471 * \param [in] minDist - minimal allowed distance from segment
472 * \param [in] curSegment - start segment
473 * \param [in] curNode - start node
474 * \param [in] segmentsOfNode - map of segments of nodes
475 * \return bool - true if a too close segment found
477 //================================================================================
479 const Segment* findTooCloseSegment( const IntPoint& p,
480 const double minDist,
482 const Segment* curSegment,
483 const SMDS_MeshNode* curNode,
484 const TSegmentsOfNode& segmentsOfNode )
486 double prevDist = Precision::Infinite();
489 double dist = SMESH_MeshAlgos::GetDistance( curSegment->myEdge, p.myNode );
490 if ( dist < minDist )
492 // check if dist is less than distance of curSegment to its cuts
493 // double minCutDist = prevDist;
494 // bool coincide = false;
495 // for ( size_t iC = 0; iC < curSegment->myCuts.size(); ++iC )
497 // if (( coincide = ( curSegment->myCuts[iC].SquareDistance( p.myNode ) < tol * tol )))
499 // for ( size_t iP = 0; iP < 2; ++iP )
501 // double cutDist = SMESH_MeshAlgos::GetDistance( curSegment->myEdge,
502 // curSegment->myCuts[iC][iP].myNode );
503 // minCutDist = std::min( minCutDist, cutDist );
506 // if ( !coincide && minCutDist > dist )
509 if ( dist > prevDist )
512 curSegment = nextSegment( curSegment, curNode, segmentsOfNode );
518 //================================================================================
520 * \brief Create a slot of given width around given 1D elements lying on a triangle mesh.
521 * The slot is constructed by cutting faces by cylindrical surfaces made around each segment.
522 * \return Edges located at the slot boundary
524 //================================================================================
526 std::vector< SMESH_MeshAlgos::Edge >
527 SMESH_MeshAlgos::MakeSlot( SMDS_ElemIteratorPtr theSegmentIt,
530 std::vector< SMDS_MeshGroup* > & theGroupsToUpdate)
532 std::vector< Edge > bndEdges;
534 if ( !theSegmentIt || !theSegmentIt->more() || !theMesh || theWidth == 0.)
537 // ----------------------------------------------------------------------------------
538 // put the input segments to a data map in order to be able finding neighboring ones
539 // ----------------------------------------------------------------------------------
541 TSegmentsOfNode segmentsOfNode;
542 ObjectPool< Segment > segmentPool;
544 while( theSegmentIt->more() )
546 const SMDS_MeshElement* edge = theSegmentIt->next();
547 if ( edge->GetType() != SMDSAbs_Edge )
548 throw SALOME_Exception( "A segment is not a mesh edge");
550 Segment* segment = segmentPool.getNew();
551 segment->myEdge = edge;
553 for ( SMDS_NodeIteratorPtr nIt = edge->nodeIterator(); nIt->more(); )
555 const SMDS_MeshNode* n = nIt->next();
556 NodeData* noData = segmentsOfNode.ChangeSeek( n );
558 noData = segmentsOfNode.Bound( n, NodeData() );
559 noData->AddSegment( segment, n );
563 // ---------------------------------
564 // Cut the mesh around the segments
565 // ---------------------------------
567 const double tol = Precision::Confusion();
568 std::vector< gp_XYZ > faceNormals;
569 SMESH_MeshAlgos::Intersector meshIntersector( theMesh, tol, faceNormals );
570 std::unique_ptr< SMESH_ElementSearcher> faceSearcher;
572 std::vector< NLink > startEdges;
573 std::vector< const SMDS_MeshNode* > faceNodes(4), edgeNodes(2);
574 std::vector<const SMDS_MeshElement *> faces(2);
575 NCollection_Map<const SMDS_MeshElement*, SMESH_Hasher > checkedFaces;
576 std::vector< IntPoint > intPoints, p(2);
577 std::vector< SMESH_NodeXYZ > facePoints(4);
578 std::vector< Intersector::TFace > cutFacePoints;
580 NCollection_DataMap< int, TGroupVec > faceID2Groups;
583 std::vector< gp_Ax1 > planeNormalVec(2);
584 gp_Ax1 * planeNormal = & planeNormalVec[0];
586 for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
588 Segment* segment = const_cast< Segment* >( segIt.next() );
590 gp_Lin segLine( segment->Ax1() );
591 gp_Ax3 cylAxis( segLine.Location(), segLine.Direction() );
592 gp_Cylinder segCylinder( cylAxis, 0.5 * theWidth );
593 double radius2( segCylinder.Radius() * segCylinder.Radius() );
595 // get normals of planes separating domains of neighboring segments
596 for ( int i = 0; i < 2; ++i ) // loop on 2 segment ends
598 const SMDS_MeshNode* n = segment->Node( i );
599 planeNormal[i] = segmentsOfNode( n ).Plane( segment );
602 // we explore faces around a segment starting from face edges;
603 // initialize a list of starting edges
606 // get a face to start searching intersected faces from
607 const SMDS_MeshNode* n0 = segment->Node( 0 );
608 SMDS_ElemIteratorPtr fIt = n0->GetInverseElementIterator( SMDSAbs_Face );
609 const SMDS_MeshElement* face = ( fIt->more() ) ? fIt->next() : 0;
610 if ( !theMesh->Contains( face ))
613 faceSearcher.reset( SMESH_MeshAlgos::GetElementSearcher( *theMesh ));
614 face = faceSearcher->FindClosestTo( SMESH_NodeXYZ( n0 ), SMDSAbs_Face );
616 // collect face edges
617 int nbNodes = face->NbCornerNodes();
618 faceNodes.assign( face->begin_nodes(), face->end_nodes() );
619 faceNodes.resize( nbNodes + 1 );
620 faceNodes[ nbNodes ] = faceNodes[ 0 ];
621 for ( int i = 0; i < nbNodes; ++i )
622 startEdges.push_back( NLink( faceNodes[i], faceNodes[i+1] ));
625 // intersect faces located around a segment
626 checkedFaces.Clear();
627 while ( !startEdges.empty() )
629 edgeNodes[0] = startEdges[0].first;
630 edgeNodes[1] = startEdges[0].second;
632 theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
633 for ( size_t iF = 0; iF < faces.size(); ++iF ) // loop on faces sharing a start edge
635 const SMDS_MeshElement* face = faces[iF];
636 if ( !checkedFaces.Add( face ))
639 int nbNodes = face->NbCornerNodes();
641 throw SALOME_Exception( "MakeSlot() accepts triangles only" );
642 faceNodes.assign( face->begin_nodes(), face->end_nodes() );
643 faceNodes.resize( nbNodes + 1 );
644 faceNodes[ nbNodes ] = faceNodes[ 0 ];
645 facePoints.assign( faceNodes.begin(), faceNodes.end() );
647 // check if cylinder axis || face
648 const gp_XYZ& faceNorm = computeNormal( face, faceNormals );
649 bool isCylinderOnFace = ( Abs( faceNorm * cylAxis.Direction().XYZ() ) < tol );
651 if ( !isCylinderOnFace )
653 if ( Intersector::CutByPlanes( face, planeNormalVec, tol, cutFacePoints ))
654 continue; // whole face cut off
655 facePoints.swap( cutFacePoints[0] );
656 facePoints.push_back( facePoints[0] );
659 // find intersection points on face edges
661 int nbPoints = facePoints.size()-1;
663 for ( int i = 0; i < nbPoints; ++i )
665 const SMESH_NodeXYZ& n1 = facePoints[i];
666 const SMESH_NodeXYZ& n2 = facePoints[i+1];
668 size_t iP = intPoints.size();
669 intersectEdge( segCylinder, n1, n2, tol, intPoints );
672 if ( isCylinderOnFace )
673 for ( ; iP < intPoints.size(); ++iP )
674 intPoints[ iP ].myEdgeIndex = i;
676 for ( ; iP < intPoints.size(); ++iP )
677 intPoints[ iP ].myEdgeIndex = edgeIndex( i, n1, n2, face,
678 intPoints[ iP ], faceNodes, tol );
680 nbFarPoints += ( segLine.SquareDistance( n1 ) > radius2 );
684 if ( nbFarPoints < nbPoints || !intPoints.empty() )
685 for ( size_t i = 1; i < faceNodes.size(); ++i )
687 const SMESH_NodeXYZ& n1 = faceNodes[i];
688 const SMESH_NodeXYZ& n2 = faceNodes[i-1];
689 isOut( n1, planeNormal, p[0].myIsOutPln );
690 isOut( n2, planeNormal, p[1].myIsOutPln );
691 if ( !isSegmentOut( p[0].myIsOutPln, p[1].myIsOutPln ))
693 startEdges.push_back( NLink( n1.Node(), n2.Node() ));
697 if ( intPoints.size() < 2 )
700 // classify intPoints by planes
701 for ( size_t i = 0; i < intPoints.size(); ++i )
702 isOut( intPoints[i].myNode, planeNormal, intPoints[i].myIsOutPln );
706 if ( intPoints.size() > 2 )
707 intPoints.push_back( intPoints[0] );
709 for ( size_t iE = 1; iE < intPoints.size(); ++iE ) // 2 <= intPoints.size() <= 5
711 if (( intPoints[iE].myIsOutPln[0] && intPoints[iE].myIsOutPln[1] ) ||
712 ( isSegmentOut( intPoints[iE].myIsOutPln, intPoints[iE-1].myIsOutPln )))
713 continue; // intPoint is out of domain
715 // check if a cutting edge connecting two intPoints is on cylinder surface
716 if ( intPoints[iE].myEdgeIndex == intPoints[iE-1].myEdgeIndex )
717 continue; // on same edge
718 if ( intPoints[iE].myNode.Node() &&
719 intPoints[iE].myNode == intPoints[iE-1].myNode ) // coincide
722 gp_XYZ edegDir = intPoints[iE].myNode - intPoints[iE-1].myNode;
724 bool toCut; // = edegDir.SquareModulus() > tol * tol;
725 if ( intPoints.size() == 2 )
727 else if ( isCylinderOnFace )
728 toCut = cylAxis.Direction().IsParallel( edegDir, tol );
731 SMESH_NodeXYZ nBetween;
732 int eInd = intPoints[iE-1].myEdgeIndex;
734 nBetween = facePoints[( 1 - (eInd-1)) % nbPoints ];
736 nBetween = faceNodes[( 1 + eInd ) % nbNodes ];
737 toCut = ( segLine.SquareDistance( nBetween ) > radius2 );
742 // limit the edge by planes
743 if ( intPoints[iE].myIsOutPln[0] ||
744 intPoints[iE].myIsOutPln[1] )
745 cutOff( intPoints[iE], intPoints[iE-1],
746 planeNormal[ intPoints[iE].myIsOutPln[1] ], tol );
748 if ( intPoints[iE-1].myIsOutPln[0] ||
749 intPoints[iE-1].myIsOutPln[1] )
750 cutOff( intPoints[iE-1], intPoints[iE],
751 planeNormal[ intPoints[iE-1].myIsOutPln[1] ], tol );
753 gp_XYZ edegDirNew = intPoints[iE].myNode - intPoints[iE-1].myNode;
754 if ( edegDir * edegDirNew < 0 ||
755 edegDir.SquareModulus() < tol * tol )
756 continue; // fully cut off
758 segment->AddCutEdge( intPoints[iE], intPoints[iE-1], face );
761 } // loop on faces sharing an edge
763 startEdges[0] = startEdges.back();
764 startEdges.pop_back();
766 } // loop on startEdges
767 } // loop on all input segments
770 // ----------------------------------------------------------
771 // If a plane fully cuts off edges of one side of a segment,
772 // it also may cut edges of adjacent segments
773 // ----------------------------------------------------------
775 for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
777 Segment* segment = const_cast< Segment* >( segIt.next() );
778 if ( segment->NbFreeEnds( tol ) >= 4 )
781 for ( int iE = 0; iE < 2; ++iE ) // loop on 2 segment ends
783 const SMDS_MeshNode* n1 = segment->Node( iE );
784 const SMDS_MeshNode* n2 = segment->Node( 1 - iE );
785 planeNormal[0] = segmentsOfNode( n1 ).Plane( segment );
788 Segment* neighborSeg = segment;
789 do // check segments connected to the segment via n2
791 neighborSeg = nextSegment( neighborSeg, n2, segmentsOfNode );
795 isNeighborCut = false;
796 for ( size_t iC = 0; iC < neighborSeg->myCuts.size(); ++iC ) // check cut edges
798 IntPoint* intPnt = &( neighborSeg->myCuts[iC].myIntPnt1 );
799 isOut( intPnt[0].myNode, planeNormal, intPnt[0].myIsOutPln, 1 );
800 isOut( intPnt[1].myNode, planeNormal, intPnt[1].myIsOutPln, 1 );
801 const Segment * closeSeg[2] = { 0, 0 };
802 if ( intPnt[0].myIsOutPln[0] )
803 closeSeg[0] = findTooCloseSegment( intPnt[0], 0.5 * theWidth - tol, tol,
804 segment, n1, segmentsOfNode );
805 if ( intPnt[1].myIsOutPln[0] )
806 closeSeg[1] = findTooCloseSegment( intPnt[1], 0.5 * theWidth - tol, tol,
807 segment, n1, segmentsOfNode );
808 int nbCut = bool( closeSeg[0] ) + bool( closeSeg[1] );
811 isNeighborCut = true;
812 if ( nbCut == 2 ) // remove a cut
814 if ( iC < neighborSeg->myCuts.size() - 1 )
815 neighborSeg->myCuts[iC] = neighborSeg->myCuts.back();
816 neighborSeg->myCuts.pop_back();
818 else // shorten cuts of 1) neighborSeg and 2) closeSeg
821 int iP = bool( closeSeg[1] );
822 gp_Lin segLine( closeSeg[iP]->Ax1() );
823 gp_Ax3 cylAxis( segLine.Location(), segLine.Direction() );
824 gp_Cylinder cyl( cylAxis, 0.5 * theWidth );
826 if ( intersectEdge( cyl, intPnt[iP].myNode, intPnt[1-iP].myNode, tol, intPoints ) &&
827 intPoints[0].SquareDistance( intPnt[iP] ) > tol * tol )
828 intPnt[iP].myNode = intPoints[0].myNode;
831 double minCutDist = theWidth;
832 gp_XYZ projection, closestProj;
834 for ( size_t iC = 0; iC < closeSeg[iP]->myCuts.size(); ++iC )
836 double cutDist = closeSeg[iP]->myCuts[iC].SquareDistance( intPnt[iP].myNode,
838 if ( cutDist < minCutDist )
840 closestProj = projection;
841 minCutDist = cutDist;
844 if ( minCutDist < tol * tol )
847 double d1 = SMESH_MeshAlgos::GetDistance( neighborSeg->myEdge,
848 closeSeg[iP]->myCuts[iCut][0].myNode );
849 double d2 = SMESH_MeshAlgos::GetDistance( neighborSeg->myEdge,
850 closeSeg[iP]->myCuts[iCut][1].myNode );
851 int iP2 = ( d2 < d1 );
852 IntPoint& ip = const_cast< IntPoint& >( closeSeg[iP]->myCuts[iCut][iP2] );
856 neighborSeg->myFreeEnds.clear();
857 neighborSeg->NbFreeEnds( tol );
860 while ( isNeighborCut );
864 // -----------------------
865 // Cut faces by cut edges
866 // -----------------------
868 for ( TSegmentIterator segIt( segmentPool ); segIt.more(); ) // loop on all segments
870 Segment* segment = const_cast< Segment* >( segIt.next() );
871 for ( size_t iC = 0; iC < segment->myCuts.size(); ++iC )
873 Cut & cut = segment->myCuts[ iC ];
874 computeNormal( cut.myFace, faceNormals );
875 meshIntersector.Cut( cut.myFace,
876 cut.myIntPnt1.myNode, cut.myIntPnt1.myEdgeIndex,
877 cut.myIntPnt2.myNode, cut.myIntPnt2.myEdgeIndex );
879 Edge e = { cut.myIntPnt1.myNode.Node(), cut.myIntPnt2.myNode.Node(), 0 };
880 bndEdges.push_back( e );
882 findGroups( cut.myFace, theGroupsToUpdate, faceID2Groups, groupVec );
886 // -----------------------------------------
887 // Make cut at the end of group of segments
888 // -----------------------------------------
890 std::vector<const SMDS_MeshElement*> polySegments;
892 for ( TSegmentsOfNode::Iterator nSegsIt( segmentsOfNode ); nSegsIt.More(); nSegsIt.Next() )
894 const NodeData& noData = nSegsIt.Value();
895 if ( noData.mySegments.size() != 1 )
898 const Segment* segment = noData.mySegments[0];
900 // find two IntPoint's of cut edges to make a cut between
901 if ( segment->myFreeEnds.size() != 4 )
902 throw SALOME_Exception( "MakeSlot(): too short end edge?" );
903 std::multimap< double, const IntPoint* > dist2IntPntMap;
904 for ( size_t iE = 0; iE < segment->myFreeEnds.size(); ++iE )
906 const SMESH_NodeXYZ& n = segment->myFreeEnds[ iE ]->myNode;
907 double d = Abs( signedDist( n, noData.myPlane ));
908 dist2IntPntMap.insert( std::make_pair( d, segment->myFreeEnds[ iE ]));
910 std::multimap< double, const IntPoint* >::iterator d2ip = dist2IntPntMap.begin();
911 SMESH_MeshAlgos::PolySegment linkNodes;
912 linkNodes.myXYZ[0] = d2ip->second->myNode;
913 linkNodes.myXYZ[1] = (++d2ip)->second->myNode;
914 linkNodes.myVector = noData.myPlane.Direction() ^ (linkNodes.myXYZ[0] - linkNodes.myXYZ[1]);
915 linkNodes.myNode1[ 0 ] = linkNodes.myNode2[ 0 ] = 0;
916 linkNodes.myNode1[ 1 ] = linkNodes.myNode2[ 1 ] = 0;
918 // create segments connecting linkNodes
919 std::vector<const SMDS_MeshElement*> newSegments;
920 std::vector<const SMDS_MeshNode*> newNodes;
921 SMESH_MeshAlgos::TListOfPolySegments polySegs(1, linkNodes);
922 SMESH_MeshAlgos::MakePolyLine( theMesh, polySegs, newSegments, newNodes,
923 /*group=*/0, faceSearcher.get() );
924 // cut faces by newSegments
926 for ( size_t i = 0; i < newSegments.size(); ++i )
928 intPoints[0].myNode = edgeNodes[0] = newSegments[i]->GetNode(0);
929 intPoints[1].myNode = edgeNodes[1] = newSegments[i]->GetNode(1);
931 // find an underlying face
932 gp_XYZ middle = 0.5 * ( intPoints[0].myNode + intPoints[1].myNode );
933 const SMDS_MeshElement* face = faceSearcher->FindClosestTo( middle, SMDSAbs_Face );
935 // find intersected edges of the face
936 int nbNodes = face->NbCornerNodes();
937 faceNodes.assign( face->begin_nodes(), face->end_nodes() );
938 faceNodes.resize( nbNodes + 1 );
939 faceNodes[ nbNodes ] = faceNodes[ 0 ];
940 for ( int iP = 0; iP < 2; ++iP )
942 intPoints[iP].myEdgeIndex = -1;
943 for ( int iN = 0; iN < nbNodes && intPoints[iP].myEdgeIndex < 0; ++iN )
945 if ( isOnEdge( faceNodes[iN], faceNodes[iN+1], intPoints[iP].myNode, tol ))
946 intPoints[iP].myEdgeIndex = iN;
952 computeNormal( face, faceNormals );
953 meshIntersector.Cut( face,
954 intPoints[0].myNode, intPoints[0].myEdgeIndex,
955 intPoints[1].myNode, intPoints[1].myEdgeIndex );
957 Edge e = { intPoints[0].myNode.Node(), intPoints[1].myNode.Node(), 0 };
958 bndEdges.push_back( e );
960 findGroups( face, theGroupsToUpdate, faceID2Groups, groupVec );
962 // add cut points to an adjacent face at ends of poly-line
963 // if they fall onto face edges
964 if (( i == 0 && intPoints[0].myEdgeIndex >= 0 ) ||
965 ( i == newSegments.size() - 1 && intPoints[1].myEdgeIndex >= 0 ))
967 for ( int iE = 0; iE < 2; ++iE ) // loop on ends of a new segment
969 if ( iE ? ( i != newSegments.size() - 1 ) : ( i != 0 ))
971 int iEdge = intPoints[ iE ].myEdgeIndex;
972 edgeNodes[0] = faceNodes[ iEdge ];
973 edgeNodes[1] = faceNodes[ iEdge+1 ];
974 theMesh->GetElementsByNodes( edgeNodes, faces, SMDSAbs_Face );
975 for ( size_t iF = 0; iF < faces.size(); ++iF )
976 if ( faces[iF] != face )
978 int iN1 = faces[iF]->GetNodeIndex( edgeNodes[0] );
979 int iN2 = faces[iF]->GetNodeIndex( edgeNodes[1] );
980 intPoints[ iE ].myEdgeIndex = Abs( iN1 - iN2 ) == 1 ? Min( iN1, iN2 ) : 2;
981 computeNormal( faces[iF], faceNormals );
982 meshIntersector.Cut( faces[iF],
983 intPoints[iE].myNode, intPoints[iE].myEdgeIndex,
984 intPoints[iE].myNode, intPoints[iE].myEdgeIndex );
986 findGroups( faces[iF], theGroupsToUpdate, faceID2Groups, groupVec );
991 } // loop on newSegments
993 polySegments.insert( polySegments.end(), newSegments.begin(), newSegments.end() );
995 } // loop on map of input segments
997 // actual mesh splitting
998 TElemIntPairVec new2OldFaces;
999 TNodeIntPairVec new2OldNodes;
1000 meshIntersector.MakeNewFaces( new2OldFaces, new2OldNodes, /*sign=*/1, /*optimize=*/true );
1002 // add new faces to theGroupsToUpdate
1003 for ( size_t i = 0; i < new2OldFaces.size(); ++i )
1005 const SMDS_MeshElement* newFace = new2OldFaces[i].first;
1006 const int oldFaceID = new2OldFaces[i].second;
1007 if ( !newFace ) continue;
1009 if ( TGroupVec* groups = const_cast< TGroupVec* >( faceID2Groups.Seek( oldFaceID )))
1010 for ( size_t iG = 0; iG < groups->size(); ++iG )
1011 (*groups)[ iG ]->Add( newFace );
1014 // remove poly-line edges
1015 for ( size_t i = 0; i < polySegments.size(); ++i )
1017 edgeNodes[0] = polySegments[i]->GetNode(0);
1018 edgeNodes[1] = polySegments[i]->GetNode(1);
1020 theMesh->RemoveFreeElement( polySegments[i] );
1022 if ( edgeNodes[0]->NbInverseElements() == 0 )
1023 theMesh->RemoveNode( edgeNodes[0] );
1024 if ( edgeNodes[1]->NbInverseElements() == 0 )
1025 theMesh->RemoveNode( edgeNodes[1] );