1 // Copyright (C) 2007-2019 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : SMESH_Offset.cxx
23 // Created : Mon Dec 25 15:52:38 2017
24 // Author : Edward AGAPOV (eap)
26 #include "SMESH_MeshAlgos.hxx"
28 #include <SMDS_PolygonalFaceOfNodes.hxx>
29 #include "SMDS_Mesh.hxx"
31 #include <Utils_SALOME_Exception.hxx>
33 #include <Bnd_B3d.hxx>
34 #include <NCollection_Map.hxx>
38 #include <boost/container/flat_set.hpp>
39 #include <boost/dynamic_bitset.hpp>
43 const int theMaxNbFaces = 256; // max number of faces sharing a node
45 typedef NCollection_DataMap< const SMDS_MeshNode*, const SMDS_MeshNode*, SMESH_Hasher > TNNMap;
46 typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap;
48 //--------------------------------------------------------------------------------
50 * \brief Intersected face side storing a node created at this intersection
51 * and a intersected face
56 const SMDS_MeshNode* myNode[2]; // side nodes. WARNING: don't set them directly, use Set()
57 mutable SMESH_NodeXYZ myIntNode; // intersection node
58 const SMDS_MeshElement* myFace; // cutter face
59 int myIndex; // index of a node on the same link
61 CutLink(const SMDS_MeshNode* node1=0,
62 const SMDS_MeshNode* node2=0,
63 const SMDS_MeshElement* face=0,
64 const int index=0) { Set ( node1, node2, face, index ); }
66 void Set( const SMDS_MeshNode* node1,
67 const SMDS_MeshNode* node2,
68 const SMDS_MeshElement* face,
71 myNode[0] = node1; myNode[1] = node2; myFace = face; myIndex = index; myReverse = false;
72 if ( myNode[0] && ( myReverse = ( myNode[0]->GetID() > myNode[1]->GetID() )))
73 std::swap( myNode[0], myNode[1] );
75 const SMDS_MeshNode* IntNode() const { return myIntNode.Node(); }
76 const SMDS_MeshNode* Node1() const { return myNode[ myReverse ]; }
77 const SMDS_MeshNode* Node2() const { return myNode[ !myReverse ]; }
79 static Standard_Integer HashCode(const CutLink& link,
80 const Standard_Integer upper)
82 Standard_Integer n = ( link.myNode[0]->GetID() +
83 link.myNode[1]->GetID() +
85 return ::HashCode( n, upper );
87 static Standard_Boolean IsEqual(const CutLink& link1, const CutLink& link2 )
89 return ( link1.myNode[0] == link2.myNode[0] &&
90 link1.myNode[1] == link2.myNode[1] &&
91 link1.myIndex == link2.myIndex );
95 typedef NCollection_Map< CutLink, CutLink > TCutLinkMap;
97 //--------------------------------------------------------------------------------
99 * \brief Part of a divided face edge
103 const SMDS_MeshNode* myNode1;
104 const SMDS_MeshNode* myNode2;
105 int myIndex; // positive -> side index, negative -> State
106 const SMDS_MeshElement* myFace;
108 enum State { _INTERNAL = -1, _COPLANAR = -2 };
110 void Set( const SMDS_MeshNode* Node1,
111 const SMDS_MeshNode* Node2,
112 const SMDS_MeshElement* Face = 0,
113 int EdgeIndex = _INTERNAL )
114 { myNode1 = Node1; myNode2 = Node2; myIndex = EdgeIndex; myFace = Face; }
116 // bool HasSameNode( const EdgePart& other ) { return ( myNode1 == other.myNode1 ||
117 // myNode1 == other.myNode2 ||
118 // myNode2 == other.myNode1 ||
119 // myNode2 == other.myNode2 );
121 bool IsInternal() const { return myIndex < 0; }
122 bool IsTwin( const EdgePart& e ) const { return myNode1 == e.myNode2 && myNode2 == e.myNode1; }
123 bool IsSame( const EdgePart& e ) const {
124 return (( myNode1 == e.myNode2 && myNode2 == e.myNode1 ) ||
125 ( myNode1 == e.myNode1 && myNode2 == e.myNode2 )); }
126 bool ReplaceCoplanar( const EdgePart& e );
127 operator SMESH_Link() const { return SMESH_Link( myNode1, myNode2 ); }
128 operator gp_Vec() const { return SMESH_NodeXYZ( myNode2 ) - SMESH_NodeXYZ( myNode1 ); }
131 //--------------------------------------------------------------------------------
133 * \brief Loop of EdgePart's forming a new face which is a part of CutFace
135 struct EdgeLoop : public SMDS_PolygonalFaceOfNodes
137 std::vector< const EdgePart* > myLinks;
138 bool myIsBndConnected; //!< is there a path to CutFace side edges
139 bool myHasPending; //!< an edge encounters twice
141 EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector<const SMDS_MeshNode *>() ) {}
142 void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; }
143 bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; }
144 bool Contains( const SMDS_MeshNode* n ) const
146 for ( size_t i = 0; i < myLinks.size(); ++i )
147 if ( myLinks[i]->myNode1 == n ) return true;
150 virtual int NbNodes() const { return myLinks.size(); }
151 virtual SMDS_ElemIteratorPtr nodesIterator() const
153 return setNodes(), SMDS_PolygonalFaceOfNodes::nodesIterator();
155 virtual SMDS_NodeIteratorPtr nodeIterator() const
157 return setNodes(), SMDS_PolygonalFaceOfNodes::nodeIterator();
159 void setNodes() const //!< set nodes to SMDS_PolygonalFaceOfNodes
161 EdgeLoop* me = const_cast<EdgeLoop*>( this );
162 me->myNodes.resize( NbNodes() );
164 for ( size_t i = 1; i < myNodes.size(); ++i ) {
165 if ( myLinks[ i ]->myNode1->GetID() < myLinks[ iMin ]->myNode1->GetID() )
168 for ( size_t i = 0; i < myNodes.size(); ++i )
169 me->myNodes[ i ] = myLinks[ ( iMin + i ) % myNodes.size() ]->myNode1;
173 //--------------------------------------------------------------------------------
175 * \brief Set of EdgeLoop's constructed from a CutFace
179 std::vector< EdgeLoop > myLoops; //!< buffer of EdgeLoop's
180 size_t myNbLoops; //!< number of constructed loops
182 const EdgePart* myEdge0; //!< & CutFace.myLinks[0]
183 size_t myNbUsedEdges; //!< nb of EdgePart's added to myLoops
184 boost::dynamic_bitset<> myIsUsedEdge; //!< is i-th EdgePart of CutFace is in any EdgeLoop
185 std::vector< EdgeLoop* > myLoopOfEdge; //!< EdgeLoop of CutFace.myLinks[i]
186 std::vector< EdgePart* > myCandidates; //!< EdgePart's starting at the same node
188 EdgeLoopSet(): myLoops(100) {}
190 void Init( const std::vector< EdgePart >& edges )
192 size_t nb = edges.size();
193 myEdge0 = & edges[0];
196 myIsUsedEdge.reset();
197 myIsUsedEdge.resize( nb, false );
198 myLoopOfEdge.clear();
199 myLoopOfEdge.resize( nb, (EdgeLoop*) 0 );
201 EdgeLoop& AddNewLoop()
203 if ( ++myNbLoops >= myLoops.size() )
204 myLoops.resize( myNbLoops + 10 );
205 myLoops[ myNbLoops-1 ].Clear();
206 return myLoops[ myNbLoops-1 ];
208 bool AllEdgesUsed() const { return myNbUsedEdges == myLoopOfEdge.size(); }
210 bool AddEdge( EdgePart& edge )
212 size_t i = Index( edge );
213 if ( myIsUsedEdge[ i ])
215 myLoops[ myNbLoops-1 ].myLinks.push_back( &edge );
216 myLoopOfEdge[ i ] = & myLoops[ myNbLoops-1 ];
217 myIsUsedEdge[ i ] = true;
221 void Erase( EdgeLoop* loop )
223 for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
224 myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0;
227 size_t Index( const EdgePart& edge ) const { return &edge - myEdge0; }
228 EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; }
231 //--------------------------------------------------------------------------------
233 * \brief Intersections of a face
237 mutable std::vector< EdgePart > myLinks;
238 const SMDS_MeshElement* myInitFace;
240 CutFace( const SMDS_MeshElement* face ): myInitFace( face ) {}
241 void AddEdge( const CutLink& p1,
243 const SMDS_MeshElement* cutter,
245 const int iNotOnPlane = -1) const;
246 void AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const;
247 bool ReplaceNodes( const TNNMap& theRm2KeepMap ) const;
249 int NbInternalEdges() const;
250 void MakeLoops( EdgeLoopSet& loops, const gp_XYZ& theFaceNorm ) const;
251 bool RemoveInternalLoops( EdgeLoopSet& theLoops ) const;
252 void CutOffLoops( EdgeLoopSet& theLoops,
253 const double theSign,
254 const std::vector< gp_XYZ >& theNormals,
255 std::vector< EdgePart >& theCutOffLinks,
256 TLinkMap& theCutOffCoplanarLinks) const;
257 void InitLinks() const;
258 bool IsCoplanar( const EdgePart* edge ) const;
260 static Standard_Integer HashCode(const CutFace& f, const Standard_Integer upper)
262 return ::HashCode( f.myInitFace->GetID(), upper );
264 static Standard_Boolean IsEqual(const CutFace& f1, const CutFace& f2 )
266 return f1.myInitFace == f2.myInitFace;
272 EdgePart* getTwin( const EdgePart* edge ) const;
275 typedef NCollection_Map< CutFace, CutFace > TCutFaceMap;
277 //--------------------------------------------------------------------------------
279 * \brief Intersection point of two edges of co-planar triangles
283 size_t myEdgeInd[2]; //!< edge indices of triangles
284 double myU [2]; //!< parameter [0,1] on edges of triangles
285 SMESH_NodeXYZ myNode; //!< intersection node
286 bool myIsCollinear;//!< edges are collinear
288 IntPoint2D() : myIsCollinear( false ) {}
290 void InitLink( CutLink& link, int iFace, const std::vector< SMESH_NodeXYZ >& nodes ) const
292 link.Set( nodes[ myEdgeInd[ iFace ] ].Node(),
293 nodes[( myEdgeInd[ iFace ] + 1 ) % nodes.size() ].Node(),
295 link.myIntNode = myNode;
297 const SMDS_MeshNode* Node() const { return myNode.Node(); }
299 struct IntPoint2DCompare
302 IntPoint2DCompare( int iFace=0 ): myI( iFace ) {}
303 bool operator() ( const IntPoint2D* ip1, const IntPoint2D* ip2 ) const
305 return ip1->myU[ myI ] < ip2->myU[ myI ];
307 bool operator() ( const IntPoint2D& ip1, const IntPoint2D& ip2 ) const
309 return ip1.myU[ myI ] < ip2.myU[ myI ];
312 typedef boost::container::flat_set< IntPoint2D, IntPoint2DCompare > TIntPointSet;
313 typedef boost::container::flat_set< IntPoint2D*, IntPoint2DCompare > TIntPointPtrSet;
315 //--------------------------------------------------------------------------------
317 * \brief Face used to find translated position of the node
321 const SMDS_MeshElement* myFace;
322 SMESH_TNodeXYZ myNode1; //!< nodes neighboring another node of myFace
323 SMESH_TNodeXYZ myNode2;
324 const gp_XYZ* myNorm;
325 bool myNodeRightOrder;
326 void operator=(const SMDS_MeshElement* f) { myFace = f; }
327 const SMDS_MeshElement* operator->() { return myFace; }
328 void SetNodes( int i0, int i1 ) //!< set myNode's
330 myNode1.Set( myFace->GetNode( i1 ));
331 int i2 = ( i0 - 1 + myFace->NbCornerNodes() ) % myFace->NbCornerNodes();
333 i2 = ( i0 + 1 ) % myFace->NbCornerNodes();
334 myNode2.Set( myFace->GetNode( i2 ));
335 myNodeRightOrder = ( Abs( i2-i1 ) == 1 ) ? i2 > i1 : i2 < i1;
337 void SetOldNodes( const SMDS_Mesh& theSrcMesh )
339 myNode1.Set( theSrcMesh.FindNode( myNode1->GetID() ));
340 myNode2.Set( theSrcMesh.FindNode( myNode2->GetID() ));
342 bool SetNormal( const std::vector< gp_XYZ >& faceNormals )
344 myNorm = & faceNormals[ myFace->GetID() ];
345 return ( myNorm->SquareModulus() > gp::Resolution() * gp::Resolution() );
347 const gp_XYZ& Norm() const { return *myNorm; }
350 //--------------------------------------------------------------------------------
352 * \brief Offset plane used to find translated position of the node
359 gp_Lin myLines[2]; //!< line of intersection with neighbor OffsetPlane's
363 void Init( const gp_XYZ& node, Face& tria, double offset )
367 myPln = gp_Pln( node + tria.Norm() * offset, tria.Norm() );
368 myIsLineOk[0] = myIsLineOk[1] = false;
369 myWeight[0] = myWeight[1] = 0;
371 bool ComputeIntersectionLine( OffsetPlane& pln );
372 void SetSkewLine( const gp_Lin& line );
373 gp_XYZ GetCommonPoint( int & nbOkPoints, double& sumWeight );
374 gp_XYZ ProjectNodeOnLine( int & nbOkPoints );
375 double Weight() const { return myWeight[0] + myWeight[1]; }
378 //================================================================================
380 * \brief Set the second line
382 //================================================================================
384 void OffsetPlane::SetSkewLine( const gp_Lin& line )
387 gp_XYZ n = myLines[0].Direction().XYZ() ^ myLines[1].Direction().XYZ();
388 if (( myIsLineOk[1] = n.SquareModulus() > gp::Resolution() ))
389 myPln = gp_Pln( myPln.Location(), n );
392 //================================================================================
394 * \brief Project myNode on myLine[0]
396 //================================================================================
398 gp_XYZ OffsetPlane::ProjectNodeOnLine( int & nbOkPoints )
400 gp_XYZ p = gp::Origin().XYZ();
403 gp_Vec l2n( myLines[0].Location(), myNode );
404 double u = l2n * myLines[0].Direction();
405 p = myLines[0].Location().XYZ() + u * myLines[0].Direction().XYZ();
411 //================================================================================
413 * \brief Computes intersection point of myLines
415 //================================================================================
417 gp_XYZ OffsetPlane::GetCommonPoint( int & nbOkPoints, double& sumWeight )
419 if ( !myIsLineOk[0] || !myIsLineOk[1] )
421 // sumWeight += myWeight[0];
422 // return ProjectNodeOnLine( nbOkPoints ) * myWeight[0];
423 return gp::Origin().XYZ();
428 gp_Vec lPerp0 = myLines[0].Direction().XYZ() ^ myPln.Axis().Direction().XYZ();
429 double dot01 = lPerp0 * myLines[1].Direction().XYZ();
430 if ( Abs( dot01 ) > 0.05 )
432 gp_Vec l0l1 = myLines[1].Location().XYZ() - myLines[0].Location().XYZ();
433 double u1 = - ( lPerp0 * l0l1 ) / dot01;
434 p = ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * u1 );
438 gp_Vec lv0( myLines[0].Location(), myNode), lv1(myLines[1].Location(), myNode );
439 double dot0( lv0 * myLines[0].Direction() ), dot1( lv1 * myLines[1].Direction() );
440 p = 0.5 * ( myLines[0].Location().XYZ() + myLines[0].Direction().XYZ() * dot0 );
441 p += 0.5 * ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * dot1 );
444 sumWeight += Weight();
450 //================================================================================
452 * \brief Compute line of intersection of 2 planes
454 //================================================================================
456 bool OffsetPlane::ComputeIntersectionLine( OffsetPlane& theNextPln )
458 const gp_XYZ& n1 = myFace->Norm();
459 const gp_XYZ& n2 = theNextPln.myFace->Norm();
461 gp_XYZ lineDir = n1 ^ n2;
464 double x = Abs( lineDir.X() );
465 double y = Abs( lineDir.Y() );
466 double z = Abs( lineDir.Z() );
468 int cooMax; // max coordinate
470 if (x > z) cooMax = 1;
474 if (y > z) cooMax = 2;
479 if ( Abs( lineDir.Coord( cooMax )) < 0.05 )
481 // parallel planes - intersection is an offset of the common edge
482 linePos = 0.5 * ( myPln.Location().XYZ() + theNextPln.myPln.Location().XYZ() );
483 lineDir = myNode - myFace->myNode2;
489 // the constants in the 2 plane equations
490 double d1 = - ( n1 * myPln.Location().XYZ() );
491 double d2 = - ( n2 * theNextPln.myPln.Location().XYZ() );
496 linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() );
497 linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() );
500 linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() );
502 linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() );
505 linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() );
506 linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() );
509 myWeight[0] = lineDir.SquareModulus();
511 myWeight[0] = 2. - myWeight[0];
513 myLines [ 0 ].SetDirection( lineDir );
514 myLines [ 0 ].SetLocation ( linePos );
515 myIsLineOk[ 0 ] = ok;
517 theNextPln.myLines [ 1 ] = myLines[ 0 ];
518 theNextPln.myIsLineOk[ 1 ] = ok;
519 theNextPln.myWeight [ 1 ] = myWeight[ 0 ];
524 //================================================================================
526 * \brief Return a translated position of a node
527 * \param [in] new2OldNodes - new and old nodes
528 * \param [in] faceNormals - normals to input faces
529 * \param [in] theSrcMesh - initial mesh
530 * \param [in] theNewPos - a computed normal
531 * \return bool - true if theNewPos is computed
533 //================================================================================
535 bool getTranslatedPosition( const SMDS_MeshNode* theNewNode,
536 const double theOffset,
538 const double theSign,
539 const std::vector< gp_XYZ >& theFaceNormals,
540 SMDS_Mesh& theSrcMesh,
543 bool useOneNormal = true;
545 // check if theNewNode needs an average position, i.e. theNewNode is convex
546 // SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
547 // const SMDS_MeshElement* f0 = faceIt->next();
548 // const gp_XYZ& norm0 = theFaceNormals[ f0->GetID() ];
549 // const SMESH_NodeXYZ nodePos = theNewNode;
550 // while ( faceIt->more() )
552 // const SMDS_MeshElement* f = faceIt->next();
553 // const int nodeInd = f->GetNodeIndex( theNewNode );
554 // SMESH_NodeXYZ nodePos2 = f->GetWrapNode( nodeInd + 1 );
556 // const gp_XYZ nnDir = ( nodePos2 - nodePos ).Normalized();
561 // const double dot = norm0 * nnDir;
566 // get faces surrounding theNewNode and sort them
567 Face faces[ theMaxNbFaces ];
568 SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator();
569 faces[0] = faceIt->next();
570 while ( !faces[0].SetNormal( theFaceNormals ) && faceIt->more() )
571 faces[0] = faceIt->next();
572 int i0 = faces[0]->GetNodeIndex( theNewNode );
573 int i1 = ( i0 + 1 ) % faces[0]->NbCornerNodes();
574 faces[0].SetNodes( i0, i1 );
575 TIDSortedElemSet elemSet, avoidSet;
577 const SMDS_MeshElement* f;
578 for ( ; faceIt->more() && iFace < theMaxNbFaces; faceIt->next() )
580 avoidSet.insert( faces[ iFace ].myFace );
581 f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
582 elemSet, avoidSet, &i0, &i1 );
585 std::reverse( &faces[0], &faces[0] + iFace + 1 );
586 for ( int i = 0; i <= iFace; ++i )
588 std::swap( faces[i].myNode1, faces[i].myNode2 );
589 faces[i].myNodeRightOrder = !faces[i].myNodeRightOrder;
591 f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(),
592 elemSet, avoidSet, &i0, &i1 );
596 faces[ ++iFace ] = f;
597 faces[ iFace ].SetNodes( i0, i1 );
598 faces[ iFace ].SetNormal( theFaceNormals );
600 int nbFaces = iFace + 1;
602 theNewPos.SetCoord( 0, 0, 0 );
603 gp_XYZ oldXYZ = SMESH_NodeXYZ( theNewNode );
605 // check if all faces are co-planar
606 bool isPlanar = true;
607 const double tol = 1e-2;
608 for ( int i = 1; i < nbFaces && isPlanar; ++i )
609 isPlanar = ( faces[i].Norm() - faces[i-1].Norm() ).SquareModulus() < tol*tol;
613 theNewPos = oldXYZ + faces[0].Norm() * theOffset;
617 // prepare OffsetPlane's
618 OffsetPlane pln[ theMaxNbFaces ];
619 for ( int i = 0; i < nbFaces; ++i )
621 faces[i].SetOldNodes( theSrcMesh );
622 pln[i].Init( oldXYZ, faces[i], theOffset );
624 // intersect neighboring OffsetPlane's
626 for ( int i = 1; i < nbFaces; ++i )
627 nbOkPoints += pln[ i-1 ].ComputeIntersectionLine( pln[ i ]);
628 nbOkPoints += pln[ nbFaces-1 ].ComputeIntersectionLine( pln[ 0 ]);
630 // move intersection lines to over parallel planes
631 if ( nbOkPoints > 1 )
632 for ( int i = 0; i < nbFaces; ++i )
633 if ( pln[i].myIsLineOk[0] && !pln[i].myIsLineOk[1] )
634 for ( int j = 1; j < nbFaces && !pln[i].myIsLineOk[1]; ++j )
636 int i2 = ( i + j ) % nbFaces;
637 if ( pln[i2].myIsLineOk[0] )
638 pln[i].SetSkewLine( pln[i2].myLines[0] );
641 // get the translated position
643 double sumWegith = 0;
644 const double minWeight = Sin( 30 * M_PI / 180. ) * Sin( 30 * M_PI / 180. );
645 for ( int i = 0; i < nbFaces; ++i )
646 if ( pln[ i ].Weight() > minWeight )
647 theNewPos += pln[ i ].GetCommonPoint( nbOkPoints, sumWegith );
649 if ( nbOkPoints == 0 )
651 // there is only one feature edge;
652 // find the theNewPos by projecting oldXYZ to any intersection line
653 for ( int i = 0; i < nbFaces; ++i )
654 theNewPos += pln[ i ].ProjectNodeOnLine( nbOkPoints );
656 if ( nbOkPoints == 0 )
658 theNewPos = oldXYZ + faces[0].Norm() * theOffset;
661 sumWegith = nbOkPoints;
663 theNewPos /= sumWegith;
666 // mark theNewNode if it is concave
667 useOneNormal = false;
668 gp_Vec moveVec( oldXYZ, theNewPos );
669 for ( int i = 0, iPrev = nbFaces-1; i < nbFaces; iPrev = i++ )
671 gp_Vec nodeVec( oldXYZ, faces[ i ].myNode1 );
672 double u = ( moveVec * nodeVec ) / nodeVec.SquareMagnitude();
673 if ( u > 0.5 ) // param [0,1] on nodeVec
675 theNewNode->setIsMarked( true );
679 gp_XYZ inFaceVec = faces[ i ].Norm() ^ nodeVec.XYZ();
680 double dot = inFaceVec * faces[ iPrev ].Norm();
681 if ( !faces[ i ].myNodeRightOrder )
683 if ( dot * theSign < 0 )
685 gp_XYZ p1 = oldXYZ + faces[ i ].Norm() * theOffset;
686 gp_XYZ p2 = oldXYZ + faces[ iPrev ].Norm() * theOffset;
687 useOneNormal = ( p1 - p2 ).SquareModulus() > 1e-12;
690 if ( useOneNormal && theNewNode->isMarked() )
699 namespace SMESH_MeshAlgos
701 //--------------------------------------------------------------------------------
703 * \brief Intersect faces of a mesh
705 struct Intersector::Algo
709 const std::vector< gp_XYZ >& myNormals;
710 TCutLinkMap myCutLinks; //!< assure sharing of new nodes
711 TCutFaceMap myCutFaces;
712 TNNMap myRemove2KeepNodes; //!< node merge map
714 // data to intersect 2 faces
715 const SMDS_MeshElement* myFace1;
716 const SMDS_MeshElement* myFace2;
717 std::vector< SMESH_NodeXYZ > myNodes1, myNodes2;
718 std::vector< double > myDist1, myDist2;
719 int myInd1, myInd2; // coordinate indices on an axis-aligned plane
720 int myNbOnPlane1, myNbOnPlane2;
721 TIntPointSet myIntPointSet;
723 Algo( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
727 //myEps( Sqrt( std::numeric_limits<double>::min() )),
728 //myEps( gp::Resolution() ),
731 void Cut( const SMDS_MeshElement* face1,
732 const SMDS_MeshElement* face2,
733 const int nbCommonNodes );
734 void Cut( const SMDS_MeshElement* face,
735 SMESH_NodeXYZ& lineEnd1,
737 SMESH_NodeXYZ& lineEnd2,
739 void MakeNewFaces( TElemIntPairVec& theNew2OldFaces,
740 TNodeIntPairVec& theNew2OldNodes,
741 const double theSign,
742 const bool theOptimize );
744 //! Cut a face by planes, whose normals point to parts to keep
745 bool CutByPlanes(const SMDS_MeshElement* face,
746 const std::vector< gp_Ax1 > & planes,
747 std::vector< SMESH_NodeXYZ > & newConnectivity );
751 bool isPlaneIntersected( const gp_XYZ& n2,
753 const std::vector< SMESH_NodeXYZ >& nodes1,
754 std::vector< double > & dist1,
757 void computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
758 const std::vector< double >& dist,
764 void addLink ( CutLink& link );
765 bool findLink( CutLink& link );
766 bool coincide( const gp_XYZ& p1, const gp_XYZ& p2, const double tol ) const
768 return ( p1 - p2 ).SquareModulus() < tol * tol;
770 gp_XY p2D( const gp_XYZ& p ) const { return gp_XY( p.Coord( myInd1 ), p.Coord( myInd2 )); }
772 void intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
773 const std::vector< double > & dist1,
775 const SMDS_MeshElement* face2,
777 void findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
778 const std::vector< double > & dist,
780 void replaceIntNode( const SMDS_MeshNode* nToKeep, const SMDS_MeshNode* nToRemove );
781 void computeIntPoint( const double u1,
786 const SMDS_MeshNode* & node1,
787 const SMDS_MeshNode* & node2);
788 void cutCollinearLink( const int iNotOnPlane1,
789 const std::vector< SMESH_NodeXYZ >& nodes1,
790 const SMDS_MeshElement* face2,
791 const CutLink& link1,
792 const CutLink& link2);
793 void setPlaneIndices( const gp_XYZ& planeNorm );
794 bool intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
795 const gp_XY s2p0, const gp_XY s2p1,
796 double & t1, double & t2,
797 bool & isCollinear );
798 bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint );
799 bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes );
800 void intersectNewEdges( const CutFace& theCFace );
801 const SMDS_MeshNode* createNode( const gp_XYZ& p );
804 //================================================================================
806 * \brief Return coordinate index with maximal abs value
808 //================================================================================
810 int MaxIndex( const gp_XYZ& x )
812 int iMaxCoo = ( Abs( x.X()) < Abs( x.Y() )) + 1;
813 if ( Abs( x.Coord( iMaxCoo )) < Abs( x.Z() ))
817 //================================================================================
819 * \brief Store a CutLink
821 //================================================================================
823 const SMDS_MeshNode* Intersector::Algo::createNode( const gp_XYZ& p )
825 const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() );
826 n->setIsMarked( true ); // cut nodes are marked
830 //================================================================================
832 * \brief Store a CutLink
834 //================================================================================
836 void Intersector::Algo::addLink( CutLink& link )
839 const CutLink* added = & myCutLinks.Added( link );
840 while ( added->myIntNode.Node() != link.myIntNode.Node() )
842 if ( !added->myIntNode )
844 added->myIntNode = link.myIntNode;
850 added = & myCutLinks.Added( link );
856 //================================================================================
858 * \brief Find a CutLink with an intersection point coincident with that of a given link
860 //================================================================================
862 bool Intersector::Algo::findLink( CutLink& link )
865 while ( myCutLinks.Contains( link ))
867 const CutLink* added = & myCutLinks.Added( link );
868 if ( !!added->myIntNode && coincide( added->myIntNode, link.myIntNode, myTol ))
870 link.myIntNode = added->myIntNode;
878 //================================================================================
880 * \brief Check if a triangle intersects the plane of another triangle
881 * \param [in] nodes1 - nodes of triangle 1
882 * \param [in] n2 - normal of triangle 2
883 * \param [in] d2 - a constant of the plane equation 2
884 * \param [out] dist1 - distance of nodes1 from the plane 2
885 * \param [out] nbOnPlane - number of nodes1 lying on the plane 2
886 * \return bool - true if the triangle intersects the plane 2
888 //================================================================================
890 bool Intersector::Algo::isPlaneIntersected( const gp_XYZ& n2,
892 const std::vector< SMESH_NodeXYZ >& nodes1,
893 std::vector< double > & dist1,
897 iNotOnPlane1 = nbOnPlane1 = 0;
898 dist1.resize( nodes1.size() );
899 for ( size_t i = 0; i < nodes1.size(); ++i )
901 dist1[i] = n2 * nodes1[i] + d2;
902 if ( Abs( dist1[i] ) < myTol )
912 if ( nbOnPlane1 == 0 )
913 for ( size_t i = 0; i < nodes1.size(); ++i )
914 if ( dist1[iNotOnPlane1] * dist1[i] < 0 )
920 //================================================================================
922 * \brief Compute parameters on the plane intersection line of intersections
923 * of edges of a triangle
924 * \param [in] nodes - triangle nodes
925 * \param [in] dist - distance of triangle nodes from the plane of another triangle
926 * \param [in] nbOnPln - number of nodes lying on the plane of another triangle
927 * \param [in] iMaxCoo - index of coordinate of max component of the plane intersection line
928 * \param [out] u - two computed parameters on the plane intersection line
929 * \param [out] iE - indices of intersected edges
931 //================================================================================
933 void Intersector::Algo::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
934 const std::vector< double >& dist,
942 u[0] = u[1] = 1e+100;
947 if ( nbOnPln == 1 && ( dist[i1] == 0. || dist[i2] == 0 ))
949 int i = dist[i1] == 0 ? i1 : i2;
950 u [ 1 ] = nodes[ i ].Coord( iMaxCoo );
954 for ( ; i2 < 3 && nb < 2; i1 = i2++ )
956 double dd = dist[i1] - dist[i2];
957 if ( dd != 0. && dist[i2] * dist[i1] <= 0. )
959 double x1 = nodes[i1].Coord( iMaxCoo );
960 double x2 = nodes[i2].Coord( iMaxCoo );
961 u [ nb ] = x1 + ( x2 - x1 ) * dist[i1] / dd;
968 std::swap( u [0], u [1] );
969 std::swap( iE[0], iE[1] );
973 //================================================================================
975 * \brief Try to find an intersection node on a link collinear with the plane intersection line
977 //================================================================================
979 void Intersector::Algo::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
980 const std::vector< double > & dist,
983 int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 );
984 CutLink link2 = link;
985 link2.Set( nodes[i1].Node(), nodes[i2].Node(), 0 );
986 if ( findLink( link2 ))
987 link.myIntNode = link2.myIntNode;
990 //================================================================================
992 * \brief Compute intersection point of a link1 with a face2
994 //================================================================================
996 void Intersector::Algo::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
997 const std::vector< double > & dist1,
999 const SMDS_MeshElement* face2,
1002 const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size();
1003 const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ];
1004 const SMESH_NodeXYZ& p2 = nodes1[ iEdge2 ];
1006 link1.Set( p1.Node(), p2.Node(), face2 );
1007 const CutLink* link = & myCutLinks.Added( link1 );
1008 if ( !link->IntNode() )
1010 if ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
1011 else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
1014 gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
1015 (gp_XYZ&)link1.myIntNode = p;
1020 gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]);
1021 while ( link->IntNode() )
1023 if ( coincide( p, link->myIntNode, myTol ))
1025 link1.myIntNode = link->myIntNode;
1029 link = & myCutLinks.Added( link1 );
1031 if ( !link1.IntNode() )
1033 if ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1;
1034 else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2;
1035 else (gp_XYZ&)link1.myIntNode = p;
1040 //================================================================================
1042 * \brief Store node replacement in myCutFaces
1044 //================================================================================
1046 void Intersector::Algo::replaceIntNode( const SMDS_MeshNode* nToKeep,
1047 const SMDS_MeshNode* nToRemove )
1049 if ( nToKeep == nToRemove )
1051 if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID
1052 myRemove2KeepNodes.Bind( nToKeep, nToRemove );
1054 myRemove2KeepNodes.Bind( nToRemove, nToKeep );
1057 //================================================================================
1059 * \brief Compute intersection point on a link of either of faces by choosing
1060 * a link whose parameter on the intersection line in maximal
1061 * \param [in] u1 - parameter on the intersection line of link iE1 of myFace1
1062 * \param [in] u2 - parameter on the intersection line of link iE2 of myFace2
1063 * \param [in] iE1 - index of a link myFace1
1064 * \param [in] iE2 - index of a link myFace2
1065 * \param [out] link - CutLink storing the intersection point
1066 * \param [out] node1 - a node of the 2nd link if two links intersect
1067 * \param [out] node2 - a node of the 2nd link if two links intersect
1069 //================================================================================
1071 void Intersector::Algo::computeIntPoint( const double u1,
1076 const SMDS_MeshNode* & node1,
1077 const SMDS_MeshNode* & node2)
1079 if ( u1 > u2 + myTol )
1081 intersectLink( myNodes1, myDist1, iE1, myFace2, link );
1083 if ( myNbOnPlane2 == 2 )
1084 findIntPointOnPlane( myNodes2, myDist2, link );
1086 else if ( u2 > u1 + myTol )
1088 intersectLink( myNodes2, myDist2, iE2, myFace1, link );
1090 if ( myNbOnPlane1 == 2 )
1091 findIntPointOnPlane( myNodes1, myDist1, link );
1093 else // edges of two faces intersect the line at the same point
1096 intersectLink( myNodes1, myDist1, iE1, myFace2, link );
1097 intersectLink( myNodes2, myDist2, iE2, myFace1, link2 );
1098 node1 = link2.Node1();
1099 node2 = link2.Node2();
1101 if ( !link.IntNode() && link2.IntNode() )
1102 link.myIntNode = link2.myIntNode;
1104 else if ( !link.IntNode() && !link2.IntNode() )
1105 (gp_XYZ&)link.myIntNode = 0.5 * ( link.myIntNode + link2.myIntNode );
1107 else if ( link.IntNode() && link2.IntNode() )
1108 replaceIntNode( link.IntNode(), link2.IntNode() );
1112 //================================================================================
1114 * \brief Add intersections to a link collinear with the intersection line
1116 //================================================================================
1118 void Intersector::Algo::cutCollinearLink( const int iNotOnPlane1,
1119 const std::vector< SMESH_NodeXYZ >& nodes1,
1120 const SMDS_MeshElement* face2,
1121 const CutLink& link1,
1122 const CutLink& link2)
1125 int iN1 = ( iNotOnPlane1 + 1 ) % 3;
1126 int iN2 = ( iNotOnPlane1 + 2 ) % 3;
1127 CutLink link( nodes1[ iN1 ].Node(), nodes1[ iN2 ].Node(), face2 );
1128 if ( link1.myFace != face2 )
1130 link.myIntNode = link1.myIntNode;
1133 if ( link2.myFace != face2 )
1135 link.myIntNode = link2.myIntNode;
1140 //================================================================================
1142 * \brief Choose indices on an axis-aligned plane
1144 //================================================================================
1146 void Intersector::Algo::setPlaneIndices( const gp_XYZ& planeNorm )
1148 switch ( MaxIndex( planeNorm )) {
1149 case 1: myInd1 = 2; myInd2 = 3; break;
1150 case 2: myInd1 = 3; myInd2 = 1; break;
1151 case 3: myInd1 = 1; myInd2 = 2; break;
1155 //================================================================================
1157 * \brief Intersect two faces
1159 //================================================================================
1161 void Intersector::Algo::Cut( const SMDS_MeshElement* face1,
1162 const SMDS_MeshElement* face2,
1163 const int nbCommonNodes)
1167 myNodes1.assign( face1->begin_nodes(), face1->end_nodes() );
1168 myNodes2.assign( face2->begin_nodes(), face2->end_nodes() );
1170 const gp_XYZ& n1 = myNormals[ face1->GetID() ];
1171 const gp_XYZ& n2 = myNormals[ face2->GetID() ];
1173 // check if triangles intersect
1174 int iNotOnPlane1, iNotOnPlane2;
1175 const double d2 = -( n2 * myNodes2[0]);
1176 if ( !isPlaneIntersected( n2, d2, myNodes1, myDist1, myNbOnPlane1, iNotOnPlane1 ))
1178 const double d1 = -( n1 * myNodes1[0]);
1179 if ( !isPlaneIntersected( n1, d1, myNodes2, myDist2, myNbOnPlane2, iNotOnPlane2 ))
1182 if ( myNbOnPlane1 == 3 || myNbOnPlane2 == 3 )// triangles are co-planar
1184 setPlaneIndices( myNbOnPlane1 == 3 ? n2 : n1 ); // choose indices on an axis-aligned plane
1187 else if ( nbCommonNodes < 2 ) // triangle planes intersect
1189 gp_XYZ lineDir = n1 ^ n2; // intersection line
1191 // check if intervals of intersections of triangles with lineDir overlap
1193 double u1[2], u2 [2]; // parameters on lineDir of edge intersection points { minU, maxU }
1194 int iE1[2], iE2[2]; // indices of edges
1195 int iMaxCoo = MaxIndex( lineDir );
1196 computeIntervals( myNodes1, myDist1, myNbOnPlane1, iMaxCoo, u1, iE1 );
1197 computeIntervals( myNodes2, myDist2, myNbOnPlane2, iMaxCoo, u2, iE2 );
1198 if ( u1[1] < u2[0] - myTol || u2[1] < u1[0] - myTol )
1199 return; // intervals do not overlap
1201 // make intersection nodes
1203 const SMDS_MeshNode *l1n1, *l1n2, *l2n1, *l2n2;
1204 CutLink link1; // intersection with smaller u on lineDir
1205 computeIntPoint( u1[0], u2[0], iE1[0], iE2[0], link1, l1n1, l1n2 );
1206 CutLink link2; // intersection with larger u on lineDir
1207 computeIntPoint( -u1[1], -u2[1], iE1[1], iE2[1], link2, l2n1, l2n2 );
1209 const CutFace& cf1 = myCutFaces.Added( CutFace( face1 ));
1210 const CutFace& cf2 = myCutFaces.Added( CutFace( face2 ));
1212 if ( coincide( link1.myIntNode, link2.myIntNode, myTol ))
1214 // intersection is a point
1215 if ( link1.IntNode() && link2.IntNode() )
1216 replaceIntNode( link1.IntNode(), link2.IntNode() );
1218 CutLink* link = link2.IntNode() ? &link2 : &link1;
1219 if ( !link->IntNode() )
1221 gp_XYZ p = 0.5 * ( link1.myIntNode + link2.myIntNode );
1222 link->myIntNode.Set( createNode( p ));
1224 if ( !link1.IntNode() ) link1.myIntNode = link2.myIntNode;
1225 if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode;
1227 cf1.AddPoint( link1, link2, myTol );
1228 cf2.AddPoint( link1, link2, myTol );
1232 // intersection is a line segment
1233 if ( !link1.IntNode() )
1234 link1.myIntNode.Set( createNode( link1.myIntNode ));
1235 if ( !link2.IntNode() )
1236 link2.myIntNode.Set( createNode( link2.myIntNode ));
1238 cf1.AddEdge( link1, link2, face2, myNbOnPlane1, iNotOnPlane1 );
1239 if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
1240 if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
1241 cf2.AddEdge( link1, link2, face1, myNbOnPlane2, iNotOnPlane2 );
1243 // add intersections to a link collinear with the intersection line
1244 if ( myNbOnPlane1 == 2 && ( link1.myFace != face2 || link2.myFace != face2 ))
1245 cutCollinearLink( iNotOnPlane1, myNodes1, face2, link1, link2 );
1247 if ( myNbOnPlane2 == 2 && ( link1.myFace != face1 || link2.myFace != face1 ))
1248 cutCollinearLink( iNotOnPlane2, myNodes2, face1, link1, link2 );
1254 } // non co-planar case
1259 //================================================================================
1261 * \brief Store a face cut by a line given by its ends
1262 * accompanied by indices of intersected face edges.
1263 * Edge index is <0 if a line end is inside the face.
1264 * \param [in] face - a face to cut
1265 * \param [inout] lineEnd1 - line end coordinates + optional node existing at this point
1266 * \param [in] edgeIndex1 - index of face edge cut by lineEnd1
1267 * \param [inout] lineEnd2 - line end coordinates + optional node existing at this point
1268 * \param [in] edgeIndex2 - index of face edge cut by lineEnd2
1270 //================================================================================
1272 void Intersector::Algo::Cut( const SMDS_MeshElement* face,
1273 SMESH_NodeXYZ& lineEnd1,
1275 SMESH_NodeXYZ& lineEnd2,
1278 if ( lineEnd1.Node() && lineEnd2.Node() &&
1279 face->GetNodeIndex( lineEnd1.Node() ) >= 0 &&
1280 face->GetNodeIndex( lineEnd2.Node() ) >= 0 )
1281 return; // intersection at a face node or edge
1283 if ((int) myNormals.size() <= face->GetID() )
1284 const_cast< std::vector< gp_XYZ >& >( myNormals ).resize( face->GetID() + 1 );
1286 const CutFace& cf = myCutFaces.Added( CutFace( face ));
1289 // look for intersection nodes coincident with line ends
1291 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1293 SMESH_NodeXYZ& lineEnd = is2nd ? lineEnd2 : lineEnd1;
1294 int edgeIndex = is2nd ? edgeIndex2 : edgeIndex1;
1295 CutLink & link = links[ is2nd ];
1297 link.myIntNode = lineEnd;
1299 for ( size_t i = ( edgeIndex < 0 ? 3 : 0 ); i < cf.myLinks.size(); ++i )
1300 if ( coincide( lineEnd, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), myTol ))
1302 link.myIntNode = cf.myLinks[i].myNode1;
1306 if ( edgeIndex >= 0 )
1308 link.Set( face->GetNode ( edgeIndex ),
1309 face->GetNodeWrap( edgeIndex + 1 ),
1314 if ( !link.myIntNode )
1315 link.myIntNode.Set( createNode( lineEnd ));
1317 lineEnd._node = link.IntNode();
1319 if ( link.myNode[0] )
1323 cf.AddEdge( links[0], links[1], /*face=*/0, /*nbOnPlane=*/0, /*iNotOnPlane=*/-1 );
1326 //================================================================================
1328 * \brief Intersect two 2D line segments
1330 //================================================================================
1332 bool Intersector::Algo::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
1333 const gp_XY s2p0, const gp_XY s2p1,
1334 double & t1, double & t2,
1335 bool & isCollinear )
1337 gp_XY u = s1p1 - s1p0;
1338 gp_XY v = s2p1 - s2p0;
1339 gp_XY w = s1p0 - s2p0;
1340 double perpDotUV = u * gp_XY( -v.Y(), v.X() );
1341 double perpDotVW = v * gp_XY( -w.Y(), w.X() );
1342 double perpDotUW = u * gp_XY( -w.Y(), w.X() );
1343 double u2 = u.SquareModulus();
1344 double v2 = v.SquareModulus();
1345 if ( u2 < myEps * myEps || v2 < myEps * myEps )
1347 if ( perpDotUV * perpDotUV / u2 / v2 < 1e-6 ) // cos ^ 2
1350 return false; // no need in collinear solution
1351 if ( perpDotUW * perpDotUW / u2 > myTol * myTol )
1352 return false; // parallel
1355 gp_XY w2 = s1p1 - s2p0;
1356 if ( Abs( v.X()) + Abs( u.X()) > Abs( v.Y()) + Abs( u.Y())) {
1357 t1 = w.X() / v.X(); // params on segment 2
1358 t2 = w2.X() / v.X();
1362 t2 = w2.Y() / v.Y();
1364 if ( Max( t1,t2 ) <= 0 || Min( t1,t2 ) >= 1 )
1365 return false; // no overlap
1368 isCollinear = false;
1370 t1 = perpDotVW / perpDotUV; // param on segment 1
1371 if ( t1 < 0. || t1 > 1. )
1372 return false; // intersection not within the segment
1374 t2 = perpDotUW / perpDotUV; // param on segment 2
1375 if ( t2 < 0. || t2 > 1. )
1376 return false; // intersection not within the segment
1381 //================================================================================
1383 * \brief Intersect two edges of co-planar triangles
1384 * \param [inout] iE1 - edge index of triangle 1
1385 * \param [inout] iE2 - edge index of triangle 2
1386 * \param [inout] intPoints - intersection points
1387 * \param [inout] nbIntPoints - nb of found intersection points
1389 //================================================================================
1391 bool Intersector::Algo::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint )
1393 int i01 = iE1, i11 = ( iE1 + 1 ) % 3;
1394 int i02 = iE2, i12 = ( iE2 + 1 ) % 3;
1395 if (( !intPoint.myIsCollinear ) &&
1396 ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
1397 myNodes1[ i01 ] == myNodes2[ i12 ] ||
1398 myNodes1[ i11 ] == myNodes2[ i02 ] ||
1399 myNodes1[ i11 ] == myNodes2[ i12 ] ))
1403 gp_XY s1p0 = p2D( myNodes1[ i01 ]);
1404 gp_XY s1p1 = p2D( myNodes1[ i11 ]);
1407 gp_XY s2p0 = p2D( myNodes2[ i02 ]);
1408 gp_XY s2p1 = p2D( myNodes2[ i12 ]);
1411 if ( !intersectEdgeEdge( s1p0,s1p1, s2p0,s2p1, t1, t2, intPoint.myIsCollinear ))
1414 intPoint.myEdgeInd[0] = iE1;
1415 intPoint.myEdgeInd[1] = iE2;
1416 intPoint.myU[0] = t1;
1417 intPoint.myU[1] = t2;
1418 (gp_XYZ&)intPoint.myNode = myNodes1[i01] * ( 1 - t1 ) + myNodes1[i11] * t1;
1420 if ( intPoint.myIsCollinear )
1423 // try to find existing node at intPoint.myNode
1425 if ( myNodes1[ i01 ] == myNodes2[ i02 ] ||
1426 myNodes1[ i01 ] == myNodes2[ i12 ] ||
1427 myNodes1[ i11 ] == myNodes2[ i02 ] ||
1428 myNodes1[ i11 ] == myNodes2[ i12 ] )
1431 const double coincTol = myTol * 1e-3;
1433 CutLink link1( myNodes1[i01].Node(), myNodes1[i11].Node(), myFace2 );
1434 CutLink link2( myNodes2[i02].Node(), myNodes2[i12].Node(), myFace1 );
1436 SMESH_NodeXYZ& n1 = myNodes1[ t1 < 0.5 ? i01 : i11 ];
1437 bool same1 = coincide( n1, intPoint.myNode, coincTol );
1440 link2.myIntNode = intPoint.myNode = n1;
1443 SMESH_NodeXYZ& n2 = myNodes2[ t2 < 0.5 ? i02 : i12 ];
1444 bool same2 = coincide( n2, intPoint.myNode, coincTol );
1447 link1.myIntNode = intPoint.myNode = n2;
1451 replaceIntNode( n1.Node(), n2.Node() );
1459 link1.myIntNode = intPoint.myNode;
1460 if ( findLink( link1 ))
1462 intPoint.myNode = link2.myIntNode = link1.myIntNode;
1467 link2.myIntNode = intPoint.myNode;
1468 if ( findLink( link2 ))
1470 intPoint.myNode = link1.myIntNode = link2.myIntNode;
1475 for ( int is2nd = 0; is2nd < 2; ++is2nd )
1477 const SMDS_MeshElement* f = is2nd ? myFace1 : myFace2;
1479 const CutFace& cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 ));
1480 for ( size_t i = 0; i < cf.myLinks.size(); ++i )
1481 if ( cf.myLinks[i].myFace == f &&
1482 //cf.myLinks[i].myIndex != EdgePart::_COPLANAR &&
1483 coincide( intPoint.myNode, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), coincTol ))
1485 intPoint.myNode.Set( cf.myLinks[i].myNode1 );
1492 intPoint.myNode._node = createNode( intPoint.myNode );
1493 link1.myIntNode = link2.myIntNode = intPoint.myNode;
1501 //================================================================================
1503 * \brief Check if a point is contained in a triangle
1505 //================================================================================
1507 bool Intersector::Algo::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes )
1510 SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ),
1511 p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ),
1513 return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
1516 //================================================================================
1518 * \brief Intersect two co-planar faces
1520 //================================================================================
1522 void Intersector::Algo::cutCoplanar()
1524 // find intersections of edges
1526 IntPoint2D intPoints[ 6 ];
1527 int nbIntPoints = 0;
1528 for ( int iE1 = 0; iE1 < 3; ++iE1 )
1530 int maxNbIntPoints = nbIntPoints + 2;
1531 for ( int iE2 = 0; iE2 < 3 && nbIntPoints < maxNbIntPoints; ++iE2 )
1532 nbIntPoints += intersectEdgeEdge( iE1, iE2, intPoints[ nbIntPoints ]);
1534 const int minNbOnPlane = Min( myNbOnPlane1, myNbOnPlane2 );
1536 if ( nbIntPoints == 0 ) // no intersections of edges
1539 if ( isPointInTriangle( myNodes1[0], myNodes2 )) // face2 includes face1
1541 else if ( isPointInTriangle( myNodes2[0], myNodes1 )) // face1 includes face2
1546 // add edges of an inner triangle to an outer one
1548 const std::vector< SMESH_NodeXYZ >& nodesIn = is1in2 ? myNodes1 : myNodes2;
1549 const SMDS_MeshElement* faceOut = is1in2 ? myFace2 : myFace1;
1550 const SMDS_MeshElement* faceIn = is1in2 ? myFace1 : myFace2;
1552 const CutFace& outFace = myCutFaces.Added( CutFace( faceOut ));
1553 CutLink link1( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
1554 CutLink link2( nodesIn.back().Node(), nodesIn.back().Node(), faceOut );
1556 link1.myIntNode = nodesIn.back();
1557 for ( size_t i = 0; i < nodesIn.size(); ++i )
1559 link2.myIntNode = nodesIn[ i ];
1560 outFace.AddEdge( link1, link2, faceIn, minNbOnPlane );
1561 link1.myIntNode = link2.myIntNode;
1566 // add parts of edges to a triangle including them
1568 CutLink link1, link2;
1569 IntPoint2D ip0, ip1;
1570 ip0.myU[0] = ip0.myU[1] = 0.;
1571 ip1.myU[0] = ip1.myU[1] = 1.;
1572 ip0.myEdgeInd[0] = ip0.myEdgeInd[1] = ip1.myEdgeInd[0] = ip1.myEdgeInd[1] = 0;
1574 for ( int isFromFace1 = 0; isFromFace1 < 2; ++isFromFace1 )
1576 const SMDS_MeshElement* faceTo = isFromFace1 ? myFace2 : myFace1;
1577 const SMDS_MeshElement* faceFrom = isFromFace1 ? myFace1 : myFace2;
1578 const std::vector< SMESH_NodeXYZ >& nodesTo = isFromFace1 ? myNodes2 : myNodes1;
1579 const std::vector< SMESH_NodeXYZ >& nodesFrom = isFromFace1 ? myNodes1 : myNodes2;
1580 const int iTo = isFromFace1 ? 1 : 0;
1581 const int iFrom = isFromFace1 ? 0 : 1;
1582 //const int nbOnPlaneFrom = isFromFace1 ? myNbOnPlane1 : myNbOnPlane2;
1584 const CutFace* cutFaceTo = & myCutFaces.Added( CutFace( faceTo ));
1585 // const CutFace* cutFaceFrom = 0;
1586 // if ( nbOnPlaneFrom > minNbOnPlane )
1587 // cutFaceFrom = & myCutFaces.Added( CutFace( faceTo ));
1589 link1.myFace = link2.myFace = faceTo;
1591 IntPoint2DCompare ipCompare( iFrom );
1592 TIntPointPtrSet pointsOnEdge( ipCompare ); // IntPoint2D sorted by parameter on edge
1594 for ( size_t iE = 0; iE < nodesFrom.size(); ++iE )
1596 // get parts of an edge iE
1598 ip0.myEdgeInd[ iTo ] = iE;
1599 ip1.myEdgeInd[ iTo ] = ( iE + 1 ) % nodesFrom.size();
1600 ip0.myNode = nodesFrom[ ip0.myEdgeInd[ iTo ]];
1601 ip1.myNode = nodesFrom[ ip1.myEdgeInd[ iTo ]];
1603 pointsOnEdge.clear();
1605 for ( int iP = 0; iP < nbIntPoints; ++iP )
1606 if ( intPoints[ iP ].myEdgeInd[ iFrom ] == iE )
1607 pointsOnEdge.insert( & intPoints[ iP ] );
1609 pointsOnEdge.insert( pointsOnEdge.begin(), & ip0 );
1610 pointsOnEdge.insert( pointsOnEdge.end(), & ip1 );
1612 // add edge parts to faceTo
1614 TIntPointPtrSet::iterator ipIt = pointsOnEdge.begin() + 1;
1615 for ( ; ipIt != pointsOnEdge.end(); ++ipIt )
1617 const IntPoint2D* p1 = *(ipIt-1);
1618 const IntPoint2D* p2 = *ipIt;
1619 gp_XYZ middle = 0.5 * ( p1->myNode + p2->myNode );
1620 if ( isPointInTriangle( middle, nodesTo ))
1622 p1->InitLink( link1, iTo, ( p1 != & ip0 ) ? nodesTo : nodesFrom );
1623 p2->InitLink( link2, iTo, ( p2 != & ip1 ) ? nodesTo : nodesFrom );
1624 cutFaceTo->AddEdge( link1, link2, faceFrom, minNbOnPlane );
1626 // if ( cutFaceFrom )
1628 // p1->InitLink( link1, iFrom, nodesFrom );
1629 // p2->InitLink( link2, iFrom, nodesFrom );
1630 // cutFaceTo->AddEdge( link1, link2, faceTo, minNbOnPlane );
1639 } // Intersector::Algo::cutCoplanar()
1641 //================================================================================
1643 * \brief Intersect edges added to myCutFaces
1645 //================================================================================
1647 void Intersector::Algo::intersectNewEdges( const CutFace& cf )
1649 IntPoint2D intPoint;
1651 if ( cf.NbInternalEdges() < 2 )
1654 if ( myNodes1.empty() )
1660 const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ];
1661 setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane
1663 size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2;
1665 for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 )
1667 if ( !cf.myLinks[i1].IsInternal() )
1670 myIntPointSet.clear();
1671 for ( size_t i2 = i1 + 2; i2 < cf.myLinks.size(); ++i2 )
1673 if ( !cf.myLinks[i2].IsInternal() )
1676 // prepare to intersection
1677 myFace1 = cf.myLinks[i1].myFace;
1678 myNodes1[0] = cf.myLinks[i1].myNode1;
1679 myNodes1[1] = cf.myLinks[i1].myNode2;
1680 myFace2 = cf.myLinks[i2].myFace;
1681 myNodes2[0] = cf.myLinks[i2].myNode1;
1682 myNodes2[1] = cf.myLinks[i2].myNode2;
1685 intPoint.myIsCollinear = true; // to find collinear solutions
1686 if ( intersectEdgeEdge( 0, 0, intPoint ))
1688 if ( cf.myLinks[i1].IsSame( cf.myLinks[i2] )) // remove i2
1690 cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i2] );
1691 cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
1695 if ( !intPoint.myIsCollinear )
1697 intPoint.myEdgeInd[1] = i2;
1698 myIntPointSet.insert( intPoint );
1700 else // if ( intPoint.myIsCollinear ) // overlapping edges
1702 myIntPointSet.clear(); // to recompute
1704 if ( intPoint.myU[0] > intPoint.myU[1] ) // orient in same direction
1706 std::swap( intPoint.myU[0], intPoint.myU[1] );
1707 std::swap( myNodes1[0], myNodes1[1] );
1709 // replace _COPLANAR by _INTERNAL
1710 cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i1+1] );
1711 cf.myLinks[i2].ReplaceCoplanar( cf.myLinks[i2+1] );
1713 if ( coincide( myNodes1[0], myNodes2[0], myTol ) &&
1714 coincide( myNodes1[1], myNodes2[1], myTol ))
1716 cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 );
1721 EdgePart common = cf.myLinks[i1];
1722 common.ReplaceCoplanar( cf.myLinks[i2] );
1724 const SMDS_MeshNode* n1 = myNodes1[0].Node(); // end nodes of an overlapping part
1725 const SMDS_MeshNode* n2 = myNodes1[1].Node();
1726 size_t i3 = cf.myLinks.size();
1728 if ( myNodes1[0] != myNodes2[0] ) // a part before the overlapping one
1730 if ( intPoint.myU[0] < 0 )
1731 cf.myLinks[i1].Set( myNodes1[0].Node(), myNodes2[0].Node(),
1732 cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
1734 cf.myLinks[i1].Set( myNodes2[0].Node(), myNodes1[0].Node(),
1735 cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
1737 cf.myLinks[i1+1].Set( cf.myLinks[i1].myNode2,
1738 cf.myLinks[i1].myNode1,
1739 cf.myLinks[i1].myFace,
1740 cf.myLinks[i1].myIndex);
1741 n1 = cf.myLinks[i1].myNode2;
1746 if ( myNodes1[1] != myNodes2[1] ) // a part after the overlapping one
1748 if ( intPoint.myU[1] < 1 )
1749 cf.myLinks[i2].Set( myNodes1[1].Node(), myNodes2[1].Node(),
1750 cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex );
1752 cf.myLinks[i2].Set( myNodes2[1].Node(), myNodes1[1].Node(),
1753 cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex );
1755 cf.myLinks[i2+1].Set( cf.myLinks[i2].myNode2,
1756 cf.myLinks[i2].myNode1,
1757 cf.myLinks[i2].myFace,
1758 cf.myLinks[i2].myIndex);
1759 n2 = cf.myLinks[i2].myNode1;
1764 if ( i3 == cf.myLinks.size() )
1765 cf.myLinks.resize( i3 + 2 );
1767 cf.myLinks[i3].Set ( n1, n2, common.myFace, common.myIndex );
1768 cf.myLinks[i3+1].Set( n2, n1, common.myFace, common.myIndex );
1770 i2 = i1 + 1; // recheck modified i1
1775 // // remember a new node
1776 // CutLink link1( myNodes1[0].Node(), myNodes1[1].Node(), cf.myInitFace );
1777 // CutLink link2( myNodes2[0].Node(), myNodes2[1].Node(), cf.myInitFace );
1778 // link2.myIntNode = link1.myIntNode = intPoint.myNode;
1779 // addLink( link1 );
1780 // addLink( link2 );
1783 // size_t i = cf.myLinks.size();
1784 // if ( intPoint.myNode != cf.myLinks[ i1 ].myNode1 &&
1785 // intPoint.myNode != cf.myLinks[ i1 ].myNode2 )
1787 // cf.myLinks.push_back( cf.myLinks[ i1 ]);
1788 // cf.myLinks.push_back( cf.myLinks[ i1 + 1 ]);
1789 // cf.myLinks[ i1 ].myNode2 = cf.myLinks[ i1 + 1 ].myNode1 = intPoint.Node();
1790 // cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = intPoint.Node();
1792 // if ( intPoint.myNode != cf.myLinks[ i2 ].myNode1 &&
1793 // intPoint.myNode != cf.myLinks[ i2 ].myNode2 )
1795 // i = cf.myLinks.size();
1796 // cf.myLinks.push_back( cf.myLinks[ i2 ]);
1797 // cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
1798 // cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = intPoint.Node();
1799 // cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = intPoint.Node();
1803 } // if ( intersectEdgeEdge( 0, 0, intPoint ))
1809 // split i1 edge and all edges it intersects
1810 // don't do it inside intersection loop in order not to loose direction of i1 edge
1811 if ( !myIntPointSet.empty() )
1813 cf.myLinks.reserve( cf.myLinks.size() + myIntPointSet.size() * 2 + 2 );
1815 EdgePart* edge1 = &cf.myLinks[ i1 ];
1816 EdgePart* twin1 = &cf.myLinks[ i1 + 1 ];
1818 TIntPointSet::iterator ipIt = myIntPointSet.begin();
1819 for ( ; ipIt != myIntPointSet.end(); ++ipIt ) // int points sorted on i1 edge
1821 size_t i = cf.myLinks.size();
1822 if ( ipIt->myNode != edge1->myNode1 &&
1823 ipIt->myNode != edge1->myNode2 )
1825 cf.myLinks.push_back( *edge1 );
1826 cf.myLinks.push_back( *twin1 );
1827 edge1->myNode2 = twin1->myNode1 = ipIt->Node();
1828 cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node();
1829 edge1 = & cf.myLinks[ i ];
1830 twin1 = & cf.myLinks[ i + 1 ];
1832 size_t i2 = ipIt->myEdgeInd[1];
1833 if ( ipIt->myNode != cf.myLinks[ i2 ].myNode1 &&
1834 ipIt->myNode != cf.myLinks[ i2 ].myNode2 )
1836 i = cf.myLinks.size();
1837 cf.myLinks.push_back( cf.myLinks[ i2 ]);
1838 cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]);
1839 cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = ipIt->Node();
1840 cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node();
1843 if ( cf.myLinks.size() >= limit )
1844 throw SALOME_Exception( "Infinite loop in Intersector::Algo::intersectNewEdges()" );
1846 ++i1; // each internal edge encounters twice
1851 //================================================================================
1853 * \brief Split intersected faces
1855 //================================================================================
1857 void Intersector::Algo::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
1858 SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
1859 const double theSign,
1860 const bool theOptimize)
1862 // fill theNew2OldFaces if empty
1863 TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin();
1864 if ( theNew2OldFaces.empty() )
1865 for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
1867 const CutFace& cf = *cutFacesIt;
1868 int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
1869 if ((int) theNew2OldFaces.size() <= index )
1870 theNew2OldFaces.resize( index + 1 );
1871 theNew2OldFaces[ index ] = std::make_pair( cf.myInitFace, index );
1874 // unmark all nodes except intersection ones
1876 for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); )
1878 const SMDS_MeshNode* n = nIt->next();
1879 if ( n->isMarked() && n->GetID()-1 < (int) theNew2OldNodes.size() )
1880 n->setIsMarked( false );
1882 // SMESH_MeshAlgos::MarkElems( myMesh->nodesIterator(), false );
1884 TCutLinkMap::const_iterator cutLinksIt = myCutLinks.cbegin();
1885 // for ( ; cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
1887 // const CutLink& link = *cutLinksIt;
1888 // if ( link.IntNode() && link.IntNode()->GetID()-1 < (int) theNew2OldNodes.size() )
1889 // link.IntNode()->setIsMarked( true );
1892 // intersect edges added to myCutFaces
1894 for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
1896 const CutFace& cf = *cutFacesIt;
1897 cf.ReplaceNodes( myRemove2KeepNodes );
1898 intersectNewEdges( cf );
1903 EdgeLoopSet loopSet;
1904 SMESH_MeshAlgos::Triangulate triangulator( theOptimize );
1905 std::vector< EdgePart > cutOffLinks;
1906 TLinkMap cutOffCoplanarLinks;
1907 std::vector< const CutFace* > touchedFaces;
1908 SMESH_MeshAlgos::TElemIntPairVec::value_type new2OldTria;
1910 std::vector< const SMDS_MeshNode* > nodes;
1911 std::vector<const SMDS_MeshElement *> faces;
1913 cutOffLinks.reserve( myCutFaces.Extent() * 2 );
1915 for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
1917 const CutFace& cf = *cutFacesIt;
1920 touchedFaces.push_back( & cf );
1924 const gp_XYZ& normal = myNormals[ cf.myInitFace->GetID() ];
1926 // form loops of new faces
1927 cf.ReplaceNodes( myRemove2KeepNodes );
1928 cf.MakeLoops( loopSet, normal );
1930 // avoid loops that are not connected to boundary edges of cf.myInitFace
1931 if ( cf.RemoveInternalLoops( loopSet ))
1933 intersectNewEdges( cf );
1934 cf.MakeLoops( loopSet, normal );
1936 // erase loops that are cut off by face intersections
1937 cf.CutOffLoops( loopSet, theSign, myNormals, cutOffLinks, cutOffCoplanarLinks );
1939 int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
1941 const SMDS_MeshElement* tria;
1942 for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
1944 EdgeLoop& loop = loopSet.myLoops[ iL ];
1945 if ( loop.myLinks.size() == 0 )
1948 int nbTria = triangulator.GetTriangles( &loop, nodes );
1949 int nbNodes = 3 * nbTria;
1950 for ( int i = 0; i < nbNodes; i += 3 )
1952 if ( nodes[i] == nodes[i+1] || nodes[i] == nodes[i+2] || nodes[i+1] == nodes[i+2] )
1955 std::cerr << "BAD tria" << std::endl;
1960 if (!( tria = myMesh->FindFace( nodes[i], nodes[i+1], nodes[i+2] )))
1961 tria = myMesh->AddFace( nodes[i], nodes[i+1], nodes[i+2] );
1962 tria->setIsMarked( true ); // not to remove it
1964 new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
1965 if ( tria->GetID() < (int)theNew2OldFaces.size() )
1966 theNew2OldFaces[ tria->GetID() ] = new2OldTria;
1968 theNew2OldFaces.push_back( new2OldTria );
1970 if ( index == tria->GetID() )
1971 index = 0; // do not remove tria
1974 theNew2OldFaces[ index ].first = 0;
1977 // remove split faces
1978 for ( size_t id = 1; id < theNew2OldFaces.size(); ++id )
1980 if ( theNew2OldFaces[id].first ||
1981 theNew2OldFaces[id].second == 0 )
1983 if ( const SMDS_MeshElement* f = myMesh->FindElement( id ))
1984 myMesh->RemoveFreeElement( f );
1987 // remove faces connected to cut off parts of cf.myInitFace
1990 for ( size_t i = 0; i < cutOffLinks.size(); ++i )
1993 nodes[0] = cutOffLinks[i].myNode1;
1994 nodes[1] = cutOffLinks[i].myNode2;
1996 if ( nodes[0] != nodes[1] &&
1997 myMesh->GetElementsByNodes( nodes, faces ))
1999 if ( cutOffLinks[i].myFace &&
2000 cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
2003 for ( size_t iF = 0; iF < faces.size(); ++iF )
2005 int index = faces[iF]->GetID();
2006 // if ( //faces[iF]->isMarked() || // kept part of cutFace
2007 // !theNew2OldFaces[ index ].first ) // already removed
2009 cutFace.myInitFace = faces[iF];
2010 // if ( myCutFaces.Contains( cutFace )) // keep cutting faces needed in CutOffLoops()
2012 // if ( !myCutFaces.Added( cutFace ).IsCut() )
2013 // theNew2OldFaces[ index ].first = 0;
2016 cutFace.myLinks.clear();
2017 cutFace.InitLinks();
2018 for ( size_t iL = 0; iL < cutFace.myLinks.size(); ++iL )
2019 if ( !cutOffLinks[i].IsSame( cutFace.myLinks[ iL ]))
2020 cutOffLinks.push_back( cutFace.myLinks[ iL ]);
2022 theNew2OldFaces[ index ].first = 0;
2023 myMesh->RemoveFreeElement( faces[iF] );
2028 // replace nodes in touched faces
2030 // treat touched faces
2031 for ( size_t i = 0; i < touchedFaces.size(); ++i )
2033 const CutFace& cf = *touchedFaces[i];
2035 int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
2036 if ( !theNew2OldFaces[ index ].first )
2037 continue; // already cut off
2039 if ( !cf.ReplaceNodes( myRemove2KeepNodes ))
2040 continue; // just keep as is
2042 if ( cf.myLinks.size() == 3 )
2044 const SMDS_MeshElement* tria = myMesh->AddFace( cf.myLinks[0].myNode1,
2045 cf.myLinks[1].myNode1,
2046 cf.myLinks[2].myNode1 );
2047 new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second );
2048 if ( tria->GetID() < (int)theNew2OldFaces.size() )
2049 theNew2OldFaces[ tria->GetID() ] = new2OldTria;
2051 theNew2OldFaces.push_back( new2OldTria );
2053 theNew2OldFaces[ index ].first = 0;
2057 // add used new nodes to theNew2OldNodes
2058 SMESH_MeshAlgos::TNodeIntPairVec::value_type new2OldNode;
2059 new2OldNode.second = 0;
2060 for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
2062 const CutLink& link = *cutLinksIt;
2063 if ( link.IntNode() ) // && link.IntNode()->NbInverseElements() > 0 )
2065 new2OldNode.first = link.IntNode();
2066 theNew2OldNodes.push_back( new2OldNode );
2073 //================================================================================
2074 Intersector::Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
2076 myAlgo = new Algo( mesh, tol, normals );
2078 //================================================================================
2079 Intersector::~Intersector()
2083 //================================================================================
2084 //! compute cut of two faces of the mesh
2085 void Intersector::Cut( const SMDS_MeshElement* face1,
2086 const SMDS_MeshElement* face2,
2087 const int nbCommonNodes )
2089 myAlgo->Cut( face1, face2, nbCommonNodes );
2091 //================================================================================
2092 //! store a face cut by a line given by its ends
2093 // accompanied by indices of intersected face edges.
2094 // Edge index is <0 if a line end is inside the face.
2095 void Intersector::Cut( const SMDS_MeshElement* face,
2096 SMESH_NodeXYZ& lineEnd1,
2098 SMESH_NodeXYZ& lineEnd2,
2101 myAlgo->Cut( face, lineEnd1, edgeIndex1, lineEnd2, edgeIndex2 );
2103 //================================================================================
2104 //! split all face intersected by Cut() methods
2105 void Intersector::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
2106 SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
2107 const double theSign,
2108 const bool theOptimize )
2110 myAlgo->MakeNewFaces( theNew2OldFaces, theNew2OldNodes, theSign, theOptimize );
2112 //================================================================================
2113 //! Cut a face by planes, whose normals point to parts to keep
2114 bool Intersector::CutByPlanes(const SMDS_MeshElement* theFace,
2115 const std::vector< gp_Ax1 > & thePlanes,
2116 const double theTol,
2117 std::vector< TFace > & theNewFaceConnectivity )
2119 theNewFaceConnectivity.clear();
2121 // check if theFace is wholly cut off
2122 std::vector< SMESH_NodeXYZ > facePoints( theFace->begin_nodes(), theFace->end_nodes() );
2123 facePoints.resize( theFace->NbCornerNodes() );
2124 for ( size_t iP = 0; iP < thePlanes.size(); ++iP )
2127 const gp_Pnt& O = thePlanes[iP].Location();
2128 for ( size_t i = 0; i < facePoints.size(); ++i )
2130 gp_Vec Op( O, facePoints[i] );
2131 nbOut += ( Op * thePlanes[iP].Direction() <= 0 );
2133 if ( nbOut == facePoints.size() )
2137 // copy theFace into a temporary mesh
2140 std::vector< const SMDS_MeshNode* > faceNodes;
2141 faceNodes.resize( facePoints.size() );
2142 for ( size_t i = 0; i < facePoints.size(); ++i )
2144 const SMESH_NodeXYZ& n = facePoints[i];
2145 faceNodes[i] = mesh.AddNode( n.X(), n.Y(), n.Z() );
2148 const SMDS_MeshElement* faceToCut = 0;
2149 switch ( theFace->NbCornerNodes() )
2152 faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2] );
2155 faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2], faceNodes[3] );
2158 faceToCut = mesh.AddPolygonalFace( faceNodes );
2161 std::vector< gp_XYZ > normals( 2 + thePlanes.size() );
2162 SMESH_MeshAlgos::FaceNormal( faceToCut, normals[ faceToCut->GetID() ]);
2164 // add faces corresponding to thePlanes
2165 std::vector< const SMDS_MeshElement* > planeFaces;
2166 double faceSize = Sqrt( faceBox.SquareExtent() );
2167 gp_XYZ center = 0.5 * ( faceBox.CornerMin() + faceBox.CornerMax() );
2168 for ( size_t i = 0; i < thePlanes.size(); ++i )
2170 gp_Ax2 plnAx( thePlanes[i].Location(), thePlanes[i].Direction() );
2171 gp_XYZ O = plnAx.Location().XYZ();
2172 gp_XYZ X = plnAx.XDirection().XYZ();
2173 gp_XYZ Y = plnAx.YDirection().XYZ();
2174 gp_XYZ Z = plnAx.Direction().XYZ();
2176 double dot = ( O - center ) * Z;
2177 gp_XYZ o = center + Z * dot; // center projected to a plane
2179 gp_XYZ p1 = o + X * faceSize * 2;
2180 gp_XYZ p2 = o + Y * faceSize * 2;
2181 gp_XYZ p3 = o - (X + Y ) * faceSize * 2;
2183 const SMDS_MeshNode* n1 = mesh.AddNode( p1.X(), p1.Y(), p1.Z() );
2184 const SMDS_MeshNode* n2 = mesh.AddNode( p2.X(), p2.Y(), p2.Z() );
2185 const SMDS_MeshNode* n3 = mesh.AddNode( p3.X(), p3.Y(), p3.Z() );
2186 planeFaces.push_back( mesh.AddFace( n1, n2, n3 ));
2188 normals[ planeFaces.back()->GetID() ] = thePlanes[i].Direction().XYZ();
2192 Algo algo ( &mesh, theTol, normals );
2193 for ( size_t i = 0; i < planeFaces.size(); ++i )
2195 algo.Cut( faceToCut, planeFaces[i], 0 );
2198 // retrieve a result
2199 SMESH_MeshAlgos::TElemIntPairVec new2OldFaces;
2200 SMESH_MeshAlgos::TNodeIntPairVec new2OldNodes;
2201 TCutFaceMap::const_iterator cutFacesIt= algo.myCutFaces.cbegin();
2202 for ( ; cutFacesIt != algo.myCutFaces.cend(); ++cutFacesIt )
2204 const CutFace& cf = *cutFacesIt;
2205 if ( cf.myInitFace != faceToCut )
2210 theNewFaceConnectivity.push_back( facePoints );
2213 // form loops of new faces
2214 EdgeLoopSet loopSet;
2215 cf.MakeLoops( loopSet, normals[ faceToCut->GetID() ]);
2217 // erase loops that are cut off by thePlanes
2218 const double sign = 1;
2219 std::vector< EdgePart > cutOffLinks;
2220 TLinkMap cutOffCoplanarLinks;
2221 cf.CutOffLoops( loopSet, sign, normals, cutOffLinks, cutOffCoplanarLinks );
2223 for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
2225 EdgeLoop& loop = loopSet.myLoops[ iL ];
2226 if ( loop.myLinks.size() > 0 )
2229 for ( SMDS_NodeIteratorPtr nIt = loop.nodeIterator(); nIt->more(); )
2231 const SMDS_MeshNode* n = nIt->next();
2232 facePoints.push_back( n );
2233 int iN = faceToCut->GetNodeIndex( n );
2235 facePoints.back()._node = 0; // an intersection point
2237 facePoints.back()._node = theFace->GetNode( iN );
2239 theNewFaceConnectivity.push_back( facePoints );
2245 return theNewFaceConnectivity.empty();
2248 } // namespace SMESH_MeshAlgos
2252 //================================================================================
2256 //================================================================================
2258 void CutFace::Dump() const
2260 std::cout << std::endl << "INI F " << myInitFace->GetID() << std::endl;
2261 for ( size_t i = 0; i < myLinks.size(); ++i )
2262 std::cout << "[" << i << "] ("
2263 << char(( myLinks[i].IsInternal() ? 'j' : '0' ) + myLinks[i].myIndex ) << ") "
2264 << myLinks[i].myNode1->GetID() << " - " << myLinks[i].myNode2->GetID()
2265 << " " << ( myLinks[i].myFace ? 'F' : 'C' )
2266 << ( myLinks[i].myFace ? myLinks[i].myFace->GetID() : 0 ) << " " << std::endl;
2269 //================================================================================
2271 * \brief Add an edge cutting this face
2272 * \param [in] p1 - start point of the edge
2273 * \param [in] p2 - end point of the edge
2274 * \param [in] cutter - a face producing the added cut edge.
2275 * \param [in] nbOnPlane - nb of triangle nodes lying on the plane of the cutter face
2277 //================================================================================
2279 void CutFace::AddEdge( const CutLink& p1,
2281 const SMDS_MeshElement* cutterFace,
2282 const int nbOnPlane,
2283 const int iNotOnPlane) const
2285 int iN[2] = { myInitFace->GetNodeIndex( p1.IntNode() ),
2286 myInitFace->GetNodeIndex( p2.IntNode() ) };
2287 if ( iN[0] >= 0 && iN[1] >= 0 )
2289 // the cutting edge is a whole side
2290 if (( cutterFace && nbOnPlane < 3 ) &&
2291 !( cutterFace->GetNodeIndex( p1.IntNode() ) >= 0 &&
2292 cutterFace->GetNodeIndex( p2.IntNode() ) >= 0 ))
2295 myLinks[ Abs( iN[0] - iN[1] ) == 1 ? Min( iN[0], iN[1] ) : 2 ].myFace = cutterFace;
2300 if ( p1.IntNode() == p2.IntNode() )
2302 AddPoint( p1, p2, 1e-10 );
2308 // cut side edges by a new one
2310 int iEOnPlane = ( nbOnPlane == 2 ) ? ( iNotOnPlane + 1 ) % 3 : -1;
2313 for ( int is2nd = 0; is2nd < 2; ++is2nd )
2315 const CutLink& p = is2nd ? p2 : p1;
2317 if ( iN[ is2nd ] >= 0 )
2320 int iE = Max( iEOnPlane, myInitFace->GetNodeIndex( p.Node1() ));
2322 continue; // link of other face
2324 SMESH_NodeXYZ n0 = myLinks[iE].myNode1;
2325 dist[ is2nd ] = ( n0 - p.myIntNode ).SquareModulus();
2327 for ( size_t i = 0; i < myLinks.size(); ++i )
2328 if ( myLinks[i].myIndex == iE )
2330 double d1 = n0.SquareDistance( myLinks[i].myNode1 );
2331 if ( d1 < dist[ is2nd ] )
2333 double d2 = n0.SquareDistance( myLinks[i].myNode2 );
2334 if ( dist[ is2nd ] < d2 )
2336 myLinks.push_back( myLinks[i] );
2337 myLinks.back().myNode1 = myLinks[i].myNode2 = p.IntNode();
2344 int state = nbOnPlane == 3 ? EdgePart::_COPLANAR : EdgePart::_INTERNAL;
2346 // look for an existing equal edge
2347 if ( nbOnPlane == 2 )
2349 SMESH_NodeXYZ n0 = myLinks[ iEOnPlane ].myNode1;
2350 if ( iN[0] >= 0 ) dist[0] = ( n0 - p1.myIntNode ).SquareModulus();
2351 if ( iN[1] >= 0 ) dist[1] = ( n0 - p2.myIntNode ).SquareModulus();
2352 if ( dist[0] > dist[1] )
2353 std::swap( dist[0], dist[1] );
2354 for ( size_t i = 0; i < myLinks.size(); ++i )
2356 if ( myLinks[i].myIndex != iEOnPlane )
2358 gp_XYZ mid = 0.5 * ( SMESH_NodeXYZ( myLinks[i].myNode1 ) +
2359 SMESH_NodeXYZ( myLinks[i].myNode2 ));
2360 double d = ( n0 - mid ).SquareModulus();
2361 if ( dist[0] < d && d < dist[1] )
2362 myLinks[i].myFace = cutterFace;
2368 EdgePart newEdge; newEdge.Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
2369 for ( size_t i = 0; i < myLinks.size(); ++i )
2371 if ( myLinks[i].IsSame( newEdge ))
2373 // if ( !myLinks[i].IsInternal() )
2374 // myLinks[ i ].myFace = cutterFace;
2376 myLinks[ i ].ReplaceCoplanar( newEdge );
2377 myLinks[ i+1 ].ReplaceCoplanar( newEdge );
2380 i += myLinks[i].IsInternal();
2384 size_t i = myLinks.size();
2385 myLinks.resize( i + 2 );
2386 myLinks[ i ].Set( p1.IntNode(), p2.IntNode(), cutterFace, state );
2387 myLinks[ i+1 ].Set( p2.IntNode(), p1.IntNode(), cutterFace, state );
2390 //================================================================================
2392 * \brief Add a point cutting this face
2394 //================================================================================
2396 void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const
2398 if ( myInitFace->GetNodeIndex( p1.IntNode() ) >= 0 ||
2399 myInitFace->GetNodeIndex( p2.IntNode() ) >= 0 )
2404 const CutLink* link = &p1;
2405 int iE = myInitFace->GetNodeIndex( link->Node1() );
2409 iE = myInitFace->GetNodeIndex( link->Node1() );
2413 // cut an existing edge by the point
2414 SMESH_NodeXYZ n0 = link->Node1();
2415 double d = ( n0 - link->myIntNode ).SquareModulus();
2417 for ( size_t i = 0; i < myLinks.size(); ++i )
2418 if ( myLinks[i].myIndex == iE )
2420 double d1 = n0.SquareDistance( myLinks[i].myNode1 );
2423 double d2 = n0.SquareDistance( myLinks[i].myNode2 );
2426 myLinks.push_back( myLinks[i] );
2427 myLinks.back().myNode1 = myLinks[i].myNode2 = link->IntNode();
2433 else // point is inside the triangle
2435 // // check if a point already added
2436 // for ( size_t i = 3; i < myLinks.size(); ++i )
2437 // if ( myLinks[i].myNode1 == p1.IntNode() )
2440 // // create a link between the point and the closest corner node
2441 // const SMDS_MeshNode* closeNode = myLinks[0].myNode1;
2442 // double minDist = p1.myIntNode.SquareDistance( closeNode );
2443 // for ( int i = 1; i < 3; ++i )
2445 // double dist = p1.myIntNode.SquareDistance( myLinks[i].myNode1 );
2446 // if ( dist < minDist )
2449 // closeNode = myLinks[i].myNode1;
2452 // if ( minDist > tol * tol )
2454 // size_t i = myLinks.size();
2455 // myLinks.resize( i + 2 );
2456 // myLinks[ i ].Set( p1.IntNode(), closeNode );
2457 // myLinks[ i+1 ].Set( closeNode, p1.IntNode() );
2462 //================================================================================
2464 * \brief Perform node replacement
2466 //================================================================================
2468 bool CutFace::ReplaceNodes( const TNNMap& theRm2KeepMap ) const
2470 bool replaced = false;
2471 for ( size_t i = 0; i < myLinks.size(); ++i )
2473 while ( theRm2KeepMap.IsBound( myLinks[i].myNode1 ))
2474 replaced = ( myLinks[i].myNode1 = theRm2KeepMap( myLinks[i].myNode1 ));
2476 while ( theRm2KeepMap.IsBound( myLinks[i].myNode2 ))
2477 replaced = ( myLinks[i].myNode2 = theRm2KeepMap( myLinks[i].myNode2 ));
2480 //if ( replaced ) // remove equal links
2482 for ( size_t i1 = 0; i1 < myLinks.size(); ++i1 )
2484 if ( myLinks[i1].myNode1 == myLinks[i1].myNode2 )
2486 myLinks.erase( myLinks.begin() + i1,
2487 myLinks.begin() + i1 + 1 + myLinks[i1].IsInternal() );
2491 size_t i2 = i1 + 1 + myLinks[i1].IsInternal();
2492 for ( ; i2 < myLinks.size(); ++i2 )
2494 if ( !myLinks[i2].IsInternal() )
2496 if ( myLinks[i1].IsSame( myLinks[i2] ))
2498 myLinks[i1]. ReplaceCoplanar( myLinks[i2] );
2499 if ( myLinks[i1].IsInternal() )
2500 myLinks[i1+1].ReplaceCoplanar( myLinks[i2+1] );
2501 if ( !myLinks[i1].myFace && myLinks[i2].myFace )
2503 myLinks[i1]. myFace = myLinks[i2].myFace;
2504 if ( myLinks[i1].IsInternal() )
2505 myLinks[i1+1].myFace = myLinks[i2+1].myFace;
2507 myLinks.erase( myLinks.begin() + i2,
2508 myLinks.begin() + i2 + 2 );
2514 i1 += myLinks[i1].IsInternal();
2521 //================================================================================
2523 * \brief Initialize myLinks with edges of myInitFace
2525 //================================================================================
2527 void CutFace::InitLinks() const
2529 if ( !myLinks.empty() ) return;
2531 int nbNodes = myInitFace->NbNodes();
2532 myLinks.reserve( nbNodes * 2 );
2533 myLinks.resize( nbNodes );
2535 for ( int i = 0; i < nbNodes; ++i )
2537 const SMDS_MeshNode* n1 = myInitFace->GetNode( i );
2538 const SMDS_MeshNode* n2 = myInitFace->GetNodeWrap( i + 1);
2539 myLinks[i].Set( n1, n2, 0, i );
2543 //================================================================================
2545 * \brief Return number of internal edges
2547 //================================================================================
2549 int CutFace::NbInternalEdges() const
2552 for ( size_t i = 3; i < myLinks.size(); ++i )
2553 nb += myLinks[i].IsInternal();
2555 return nb / 2; // each internal edge encounters twice
2558 //================================================================================
2560 * \brief Remove loops that are not connected to boundary edges of myFace by
2561 * adding edges connecting these loops to the boundary
2563 //================================================================================
2565 bool CutFace::RemoveInternalLoops( EdgeLoopSet& theLoops ) const
2567 size_t nbReachedLoops = 0;
2569 // count loops including boundary EdgeParts
2570 for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2572 EdgeLoop& loop = theLoops.myLoops[ iL ];
2574 for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2575 if ( !loop.myLinks[ iE ]->IsInternal() )
2577 nbReachedLoops += loop.SetConnected();
2581 if ( nbReachedLoops == theLoops.myNbLoops )
2582 return false; // no unreachable loops
2585 // try to reach all loops by propagating via internal edges shared by loops
2586 size_t prevNbReached;
2589 prevNbReached = nbReachedLoops;
2591 for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2593 EdgeLoop& loop = theLoops.myLoops[ iL ];
2594 if ( !loop.myIsBndConnected )
2597 for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2598 if ( loop.myLinks[ iE ]->IsInternal() )
2600 const EdgePart* twinEdge = getTwin( loop.myLinks[ iE ]);
2601 EdgeLoop* loop2 = theLoops.GetLoopOf( twinEdge );
2602 if ( loop2->SetConnected() && ++nbReachedLoops == theLoops.myNbLoops )
2603 return false; // no unreachable loops
2607 while ( prevNbReached < nbReachedLoops );
2610 // add links connecting internal loops with the boundary ones
2612 for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
2614 EdgeLoop& loop = theLoops.myLoops[ iL ];
2615 if ( loop.myIsBndConnected )
2618 // find a pair of closest nodes
2619 const SMDS_MeshNode *closestNode1, *closestNode2;
2620 double minDist = 1e100;
2621 for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE )
2623 SMESH_NodeXYZ n1 = loop.myLinks[ iE ]->myNode1;
2625 for ( size_t i = 0; i < myLinks.size(); ++i )
2627 if ( !loop.Contains( myLinks[i].myNode1 ))
2629 double dist = n1.SquareDistance( myLinks[i].myNode1 );
2630 if ( dist < minDist )
2633 closestNode1 = loop.myLinks[ iE ]->myNode1;
2634 closestNode2 = myLinks[i].myNode1;
2637 if ( myLinks[i].IsInternal() )
2642 size_t i = myLinks.size();
2643 myLinks.resize( i + 2 );
2644 myLinks[ i ].Set( closestNode1, closestNode2 );
2645 myLinks[ i+1 ].Set( closestNode2, closestNode1 );
2651 //================================================================================
2653 * \brief Return equal reversed edge
2655 //================================================================================
2657 EdgePart* CutFace::getTwin( const EdgePart* edge ) const
2659 size_t i = edge - & myLinks[0];
2661 if ( i > 2 && myLinks[ i-1 ].IsTwin( *edge ))
2662 return & myLinks[ i-1 ];
2664 if ( i+1 < myLinks.size() &&
2665 myLinks[ i+1 ].IsTwin( *edge ))
2666 return & myLinks[ i+1 ];
2671 //================================================================================
2673 * \brief Fill loops of edges
2675 //================================================================================
2677 void CutFace::MakeLoops( EdgeLoopSet& theLoops, const gp_XYZ& theFaceNorm ) const
2679 theLoops.Init( myLinks );
2681 if ( myLinks.size() == 3 )
2683 theLoops.AddNewLoop();
2684 theLoops.AddEdge( myLinks[0] );
2685 theLoops.AddEdge( myLinks[1] );
2686 theLoops.AddEdge( myLinks[2] );
2690 while ( !theLoops.AllEdgesUsed() )
2692 theLoops.AddNewLoop();
2694 // add 1st edge to a new loop
2696 for ( i1 = theLoops.myNbLoops - 1; i1 < myLinks.size(); ++i1 )
2697 if ( theLoops.AddEdge( myLinks[i1] ))
2700 EdgePart* lastEdge = & myLinks[ i1 ];
2701 EdgePart* twinEdge = getTwin( lastEdge );
2702 const SMDS_MeshNode* firstNode = lastEdge->myNode1;
2703 const SMDS_MeshNode* lastNode = lastEdge->myNode2;
2705 do // add the rest edges
2707 theLoops.myCandidates.clear(); // edges starting at lastNode
2710 // find candidate edges
2711 for ( size_t i = i1 + 1; i < myLinks.size(); ++i )
2712 if ( myLinks[ i ].myNode1 == lastNode &&
2713 &myLinks[ i ] != twinEdge &&
2714 !theLoops.myIsUsedEdge[ i ])
2716 theLoops.myCandidates.push_back( & myLinks[ i ]);
2717 nbInternal += myLinks[ i ].IsInternal();
2720 // choose among candidates
2721 if ( theLoops.myCandidates.size() == 0 )
2723 theLoops.GetLoopOf( lastEdge )->myHasPending = true;
2724 lastEdge = twinEdge;
2726 else if ( theLoops.myCandidates.size() == 1 )
2728 lastEdge = theLoops.myCandidates[0];
2730 else if ( nbInternal == 1 && !lastEdge->IsInternal() )
2732 lastEdge = theLoops.myCandidates[ !theLoops.myCandidates[0]->IsInternal() ];
2736 gp_Vec lastVec = *lastEdge;
2737 double maxAngle = -2 * M_PI;
2738 for ( size_t i = 0; i < theLoops.myCandidates.size(); ++i )
2740 double angle = lastVec.AngleWithRef( *theLoops.myCandidates[i], theFaceNorm );
2741 if ( angle > maxAngle )
2744 lastEdge = theLoops.myCandidates[i];
2748 theLoops.AddEdge( *lastEdge );
2749 lastNode = lastEdge->myNode2;
2750 twinEdge = getTwin( lastEdge );
2752 while ( lastNode != firstNode );
2754 } // while ( !theLoops.AllEdgesUsed() )
2759 //================================================================================
2761 * \brief Erase loops that are cut off by face intersections
2763 //================================================================================
2765 void CutFace::CutOffLoops( EdgeLoopSet& theLoops,
2766 const double theSign,
2767 const std::vector< gp_XYZ >& theNormals,
2768 std::vector< EdgePart >& theCutOffLinks,
2769 TLinkMap& theCutOffCoplanarLinks) const
2772 for ( size_t i = 0; i < myLinks.size(); ++i )
2774 if ( !myLinks[i].myFace )
2777 EdgeLoop* loop = theLoops.GetLoopOf( & myLinks[i] );
2778 if ( !loop || loop->myLinks.empty() || loop->myHasPending )
2781 bool toErase, isCoplanar = ( myLinks[i].myIndex == EdgePart::_COPLANAR );
2783 gp_Vec iniNorm = theNormals[ myInitFace->GetID() ];
2786 toErase = ( myLinks[i].myFace->GetID() > myInitFace->GetID() );
2788 const EdgePart* twin = getTwin( & myLinks[i] );
2789 if ( !twin || twin->myFace == myLinks[i].myFace )
2791 // only one co-planar face includes myLinks[i]
2792 gp_Vec inFaceDir = iniNorm ^ myLinks[i];
2793 gp_XYZ edgePnt = SMESH_NodeXYZ( myLinks[i].myNode1 );
2794 for ( int iN = 0; iN < 3; ++iN )
2796 gp_Vec inCutFaceDir = ( SMESH_NodeXYZ( myLinks[i].myFace->GetNode( iN )) - edgePnt );
2797 if ( inCutFaceDir * inFaceDir < 0 )
2807 gp_Vec cutNorm = theNormals[ myLinks[i].myFace->GetID() ];
2808 gp_Vec inFaceDir = iniNorm ^ myLinks[i];
2810 toErase = inFaceDir * cutNorm * theSign < 0;
2813 // erase a neighboring loop
2815 if ( const EdgePart* twin = getTwin( & myLinks[i] ))
2816 loop = theLoops.GetLoopOf( twin );
2817 toErase = ( loop && !loop->myLinks.empty() );
2825 // remember whole sides of myInitFace that are cut off
2826 for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE )
2828 if ( !loop->myLinks[ iE ]->myFace &&
2829 !loop->myLinks[ iE ]->IsInternal() )// &&
2830 // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
2831 // !loop->myLinks[ iE ]->myNode2->isMarked() )
2833 int i = loop->myLinks[ iE ]->myIndex;
2834 sideEdge.Set( myInitFace->GetNode ( i ),
2835 myInitFace->GetNodeWrap( i+1 ));
2836 theCutOffLinks.push_back( sideEdge );
2838 if ( !sideEdge.IsSame( *loop->myLinks[ iE ] )) // nodes replaced
2840 theCutOffLinks.push_back( *loop->myLinks[ iE ] );
2843 else if ( IsCoplanar( loop->myLinks[ iE ]))
2845 // propagate erasure to a co-planar face
2846 theCutOffLinks.push_back( *loop->myLinks[ iE ]);
2848 else if ( loop->myLinks[ iE ]->myFace &&
2849 loop->myLinks[ iE ]->IsInternal() )
2850 theCutOffLinks.push_back( *loop->myLinks[ iE ]);
2854 theLoops.Erase( loop );
2861 //================================================================================
2863 * \brief Check if the face has cut edges
2865 //================================================================================
2867 bool CutFace::IsCut() const
2869 if ( myLinks.size() > 3 )
2872 if ( myLinks.size() == 3 )
2873 for ( size_t i = 0; i < 3; ++i )
2874 if ( myLinks[i].myFace )
2880 //================================================================================
2882 * \brief Check if an edge is produced by a co-planar cut
2884 //================================================================================
2886 bool CutFace::IsCoplanar( const EdgePart* edge ) const
2888 if ( edge->myIndex == EdgePart::_COPLANAR )
2890 const EdgePart* twin = getTwin( edge );
2891 return ( !twin || twin->myIndex == EdgePart::_COPLANAR );
2896 //================================================================================
2898 * \brief Replace _COPLANAR cut edge by _INTERNAL or vice versa
2900 //================================================================================
2902 bool EdgePart::ReplaceCoplanar( const EdgePart& e )
2904 if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL )
2906 //check if the faces are connected
2907 int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size();
2908 bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
2909 ( myIndex == _COPLANAR && nbCommonNodes < 2 ));
2912 myIndex = e.myIndex;
2922 //================================================================================
2924 * \brief Create an offsetMesh of given faces
2925 * \param [in] faceIt - the input faces
2926 * \param [out] new2OldFaces - history of faces (new face -> old face ID)
2927 * \param [out] new2OldNodes - history of nodes (new node -> old node ID)
2928 * \return SMDS_Mesh* - the new offset mesh, a caller should delete
2930 //================================================================================
2932 SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt,
2933 SMDS_Mesh& theSrcMesh,
2934 const double theOffset,
2935 const bool theFixIntersections,
2936 TElemIntPairVec& theNew2OldFaces,
2937 TNodeIntPairVec& theNew2OldNodes)
2939 if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 )
2940 throw SALOME_Exception( "Offset of quadratic mesh not supported" );
2941 if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() )
2942 throw SALOME_Exception( "Offset of non-triangular mesh not supported" );
2944 SMDS_Mesh* newMesh = new SMDS_Mesh;
2945 theNew2OldFaces.clear();
2946 theNew2OldNodes.clear();
2947 theNew2OldFaces.push_back
2948 ( std::make_pair(( const SMDS_MeshElement*) 0, 0)); // to have index == face->GetID()
2950 // copy input faces to the newMesh keeping IDs of nodes
2952 double minNodeDist = 1e100;
2954 std::vector< const SMDS_MeshNode* > nodes;
2955 while ( theFaceIt->more() )
2957 const SMDS_MeshElement* face = theFaceIt->next();
2958 if ( face->GetType() != SMDSAbs_Face ) continue;
2961 nodes.assign( face->begin_nodes(), face->end_nodes() );
2962 for ( size_t i = 0; i < nodes.size(); ++i )
2964 const SMDS_MeshNode* newNode = newMesh->FindNode( nodes[i]->GetID() );
2967 SMESH_NodeXYZ xyz( nodes[i] );
2968 newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() );
2969 theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i]->GetID() ));
2973 const SMDS_MeshElement* newFace = 0;
2974 switch ( face->GetEntityType() )
2976 case SMDSEntity_Triangle:
2977 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2] );
2979 case SMDSEntity_Quad_Triangle:
2980 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
2981 nodes[3],nodes[4],nodes[5] );
2983 case SMDSEntity_BiQuad_Triangle:
2984 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],
2985 nodes[3],nodes[4],nodes[5],nodes[6] );
2987 case SMDSEntity_Quadrangle:
2988 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3] );
2990 case SMDSEntity_Quad_Quadrangle:
2991 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],
2992 nodes[4],nodes[5],nodes[6],nodes[7] );
2994 case SMDSEntity_BiQuad_Quadrangle:
2995 newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4],
2996 nodes[5],nodes[6],nodes[7],nodes[8] );
2998 case SMDSEntity_Polygon:
2999 newFace = newMesh->AddPolygonalFace( nodes );
3001 case SMDSEntity_Quad_Polygon:
3002 newFace = newMesh->AddQuadPolygonalFace( nodes );
3007 theNew2OldFaces.push_back( std::make_pair( newFace, face->GetID() ));
3009 SMESH_NodeXYZ pPrev = nodes.back(), p;
3010 for ( size_t i = 0; i < nodes.size(); ++i )
3013 double dist = ( pPrev - p ).SquareModulus();
3014 if ( dist > std::numeric_limits<double>::min() )
3018 } // while ( faceIt->more() )
3021 // compute normals to faces
3022 std::vector< gp_XYZ > normals( theNew2OldFaces.size() );
3023 for ( size_t i = 1; i < normals.size(); ++i )
3025 if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].first, normals[i] ))
3026 normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
3029 const double sign = ( theOffset < 0 ? -1 : +1 );
3030 const double tol = Min( 1e-3 * Sqrt( minNodeDist ),
3031 1e-2 * theOffset * sign );
3033 // translate new nodes by normal to input faces
3035 std::vector< const SMDS_MeshNode* > multiNormalNodes;
3036 for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
3038 const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
3040 if ( getTranslatedPosition( newNode, theOffset, tol*10., sign, normals, theSrcMesh, newXYZ ))
3041 newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3043 multiNormalNodes.push_back( newNode );
3045 // make multi-normal translation
3046 std::vector< SMESH_NodeXYZ > multiPos(10);
3047 for ( size_t i = 0; i < multiNormalNodes.size(); ++i )
3049 const SMDS_MeshNode* newNode = multiNormalNodes[i];
3050 newNode->setIsMarked( true );
3051 SMESH_NodeXYZ oldXYZ = newNode;
3053 for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3055 const SMDS_MeshElement* newFace = fIt->next();
3056 const int faceIndex = newFace->GetID();
3057 const gp_XYZ& oldNorm = normals[ faceIndex ];
3058 const gp_XYZ newXYZ = oldXYZ + oldNorm * theOffset;
3059 if ( multiPos.empty() )
3061 newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3062 multiPos.emplace_back( newNode );
3067 for ( size_t iP = 0; iP < multiPos.size() && !newNode; ++iP )
3068 if (( multiPos[iP] - newXYZ ).SquareModulus() < tol * tol )
3069 newNode = multiPos[iP].Node();
3072 newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
3073 newNode->setIsMarked( true );
3074 theNew2OldNodes.push_back( std::make_pair( newNode, 0 ));
3075 multiPos.emplace_back( newNode );
3078 if ( newNode != oldXYZ.Node() )
3080 nodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
3081 nodes[ newFace->GetNodeIndex( oldXYZ.Node() )] = newNode;
3082 newMesh->ChangeElementNodes( newFace, & nodes[0], nodes.size() );
3087 if ( !theFixIntersections )
3091 // remove new faces around concave nodes (they are marked) if the faces are inverted
3093 for ( size_t i = 0; i < theNew2OldNodes.size(); ++i )
3095 const SMDS_MeshNode* newNode = theNew2OldNodes[i].first;
3096 //const SMDS_MeshNode* oldNode = theNew2OldNodes[i].second;
3097 if ( newNode->isMarked() )
3099 //gp_XYZ moveVec = sign * ( SMESH_NodeXYZ( newNode ) - SMESH_NodeXYZ( oldNode ));
3101 //bool haveInverseFace = false;
3102 for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3104 const SMDS_MeshElement* newFace = fIt->next();
3105 const int faceIndex = newFace->GetID();
3106 const gp_XYZ& oldNorm = normals[ faceIndex ];
3107 if ( !SMESH_MeshAlgos::FaceNormal( newFace, faceNorm, /*normalize=*/false ) ||
3108 //faceNorm * moveVec < 0 )
3109 faceNorm * oldNorm < 0 )
3111 //haveInverseFace = true;
3112 theNew2OldFaces[ faceIndex ].first = 0;
3113 newMesh->RemoveFreeElement( newFace );
3117 // if ( haveInverseFace )
3119 // newMesh->MoveNode( newNode, oldNode->X(), oldNode->Y(), oldNode->Z() );
3121 // for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); )
3123 // const SMDS_MeshElement* newFace = fIt->next();
3124 // if ( !SMESH_MeshAlgos::FaceNormal( newFace, normals[ newFace->GetID() ] ))
3125 // normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
3129 // mark all new nodes located closer than theOffset from theSrcMesh
3132 // ==================================================
3133 // find self-intersections of new faces and fix them
3134 // ==================================================
3136 std::unique_ptr< SMESH_ElementSearcher > fSearcher
3137 ( SMESH_MeshAlgos::GetElementSearcher( *newMesh, tol ));
3139 Intersector intersector( newMesh, tol, normals );
3141 std::vector< const SMDS_MeshElement* > closeFaces;
3142 std::vector< const SMDS_MeshNode* > faceNodes;
3144 for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF )
3146 const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first;
3147 if ( !newFace ) continue;
3148 faceNodes.assign( newFace->begin_nodes(), newFace->end_nodes() );
3150 bool isConcaveNode1 = false;
3151 for ( size_t iN = 0; iN < faceNodes.size() && !isConcaveNode1; ++iN )
3152 isConcaveNode1 = faceNodes[iN]->isMarked();
3154 // get faces close to a newFace
3157 for ( size_t i = 0; i < faceNodes.size(); ++i )
3158 faceBox.Add( SMESH_NodeXYZ( faceNodes[i] ));
3159 faceBox.Enlarge( tol );
3161 fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces );
3163 // intersect the newFace with closeFaces
3165 for ( size_t i = 0; i < closeFaces.size(); ++i )
3167 const SMDS_MeshElement* closeFace = closeFaces[i];
3168 if ( closeFace->GetID() <= newFace->GetID() )
3169 continue; // this pair already treated
3171 // do not intersect connected faces if they have no concave nodes
3172 int nbCommonNodes = 0;
3173 for ( size_t iN = 0; iN < faceNodes.size(); ++iN )
3174 nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN] ) >= 0 );
3176 if ( !isConcaveNode1 )
3178 bool isConcaveNode2 = false;
3179 for ( SMDS_ElemIteratorPtr nIt = closeFace->nodesIterator(); nIt->more(); )
3180 if (( isConcaveNode2 = nIt->next()->isMarked() ))
3183 if ( !isConcaveNode2 && nbCommonNodes > 0 )
3187 intersector.Cut( newFace, closeFace, nbCommonNodes );
3190 intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign );