-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2020 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
{
const int theMaxNbFaces = 256; // max number of faces sharing a node
- typedef NCollection_DataMap< Standard_Address, const SMDS_MeshNode* > TNNMap;
- typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap;
+ typedef NCollection_DataMap< const SMDS_MeshNode*, const SMDS_MeshNode*, SMESH_Hasher > TNNMap;
+ typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap;
//--------------------------------------------------------------------------------
/*!
* \brief Intersected face side storing a node created at this intersection
- * and a intersected face
+ * and an intersected face
*/
struct CutLink
{
bool myReverse;
- const SMDS_MeshNode* myNode[2]; // side nodes
+ const SMDS_MeshNode* myNode[2]; // side nodes. WARNING: don't set them directly, use Set()
mutable SMESH_NodeXYZ myIntNode; // intersection node
const SMDS_MeshElement* myFace; // cutter face
int myIndex; // index of a node on the same link
int myIndex; // positive -> side index, negative -> State
const SMDS_MeshElement* myFace;
- enum State { _INTERNAL = -1, _COPLANAR = -2 };
+ enum State { _INTERNAL = -1, _COPLANAR = -2, _PENDING = -3 };
void Set( const SMDS_MeshNode* Node1,
const SMDS_MeshNode* Node2,
EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector<const SMDS_MeshNode *>() ) {}
void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; }
bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; }
- bool Contains( const SMDS_MeshNode* n ) const
+ size_t Contains( const SMDS_MeshNode* n ) const
{
for ( size_t i = 0; i < myLinks.size(); ++i )
- if ( myLinks[i]->myNode1 == n ) return true;
- return false;
+ if ( myLinks[i]->myNode1 == n ) return i + 1;
+ return 0;
}
virtual int NbNodes() const { return myLinks.size(); }
virtual SMDS_ElemIteratorPtr nodesIterator() const
myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0;
loop->Clear();
}
+ void Join( EdgeLoop& loop1, size_t iAfterConcact,
+ EdgeLoop& loop2, size_t iFromEdge2 )
+ {
+ std::vector< const EdgePart* > linksAfterContact( loop1.myLinks.begin() + iAfterConcact,
+ loop1.myLinks.end() );
+ loop1.myLinks.reserve( loop2.myLinks.size() + loop1.myLinks.size() );
+ loop1.myLinks.resize( iAfterConcact );
+ loop1.myLinks.insert( loop1.myLinks.end(),
+ loop2.myLinks.begin() + iFromEdge2, loop2.myLinks.end() );
+ loop1.myLinks.insert( loop1.myLinks.end(),
+ loop2.myLinks.begin(), loop2.myLinks.begin() + iFromEdge2 );
+ loop1.myLinks.insert( loop1.myLinks.end(),
+ linksAfterContact.begin(), linksAfterContact.end() );
+ loop1.myIsBndConnected = loop2.myIsBndConnected;
+ loop2.Clear();
+ for ( size_t iE = 0; iE < loop1.myLinks.size(); ++iE )
+ myLoopOfEdge[ Index( *loop1.myLinks[ iE ] )] = & loop1;
+ }
size_t Index( const EdgePart& edge ) const { return &edge - myEdge0; }
EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; }
};
return useOneNormal;
}
+ //================================================================================
+ /*!
+ * \brief Remove small faces
+ */
+ //================================================================================
+
+ void removeSmallFaces( SMDS_Mesh* theMesh,
+ SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
+ const double theTol2 )
+ {
+ std::vector< SMESH_NodeXYZ > points(3);
+ std::vector< const SMDS_MeshNode* > nodes(3);
+ for ( SMDS_ElemIteratorPtr faceIt = theMesh->elementsIterator(); faceIt->more(); )
+ {
+ const SMDS_MeshElement* face = faceIt->next();
+ points.assign( face->begin_nodes(), face->end_nodes() );
+
+ SMESH_NodeXYZ* prevN = & points.back();
+ for ( size_t i = 0; i < points.size(); ++i )
+ {
+ double dist2 = ( *prevN - points[ i ]).SquareModulus();
+ if ( dist2 < theTol2 )
+ {
+ const SMDS_MeshNode* nToRemove =
+ (*prevN)->GetID() > points[ i ]->GetID() ? prevN->Node() : points[ i ].Node();
+ const SMDS_MeshNode* nToKeep =
+ nToRemove == points[ i ].Node() ? prevN->Node() : points[ i ].Node();
+ for ( SMDS_ElemIteratorPtr fIt = nToRemove->GetInverseElementIterator(); fIt->more(); )
+ {
+ const SMDS_MeshElement* f = fIt->next();
+ if ( f == face )
+ continue;
+ nodes.assign( f->begin_nodes(), f->end_nodes() );
+ nodes[ f->GetNodeIndex( nToRemove )] = nToKeep;
+ theMesh->ChangeElementNodes( f, &nodes[0], nodes.size() );
+ }
+ theNew2OldFaces[ face->GetID() ].first = 0;
+ theMesh->RemoveFreeElement( face );
+ break;
+ }
+ prevN = & points[ i ];
+ }
+ continue;
+ }
+ return;
+ }
+
+} // namespace
+
+namespace SMESH_MeshAlgos
+{
//--------------------------------------------------------------------------------
/*!
* \brief Intersect faces of a mesh
*/
- struct Intersector
+ struct Intersector::Algo
{
SMDS_Mesh* myMesh;
double myTol, myEps;
int myNbOnPlane1, myNbOnPlane2;
TIntPointSet myIntPointSet;
- Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
+ Algo( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
: myMesh( mesh ),
myTol( tol ),
myEps( 1e-100 ),
void Cut( const SMDS_MeshElement* face1,
const SMDS_MeshElement* face2,
const int nbCommonNodes );
- void MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces,
- SMESH_MeshAlgos::TNPairVec& theNew2OldNodes,
- const double theSign );
+ void Cut( const SMDS_MeshElement* face,
+ SMESH_NodeXYZ& lineEnd1,
+ int edgeIndex1,
+ SMESH_NodeXYZ& lineEnd2,
+ int edgeIndex2 );
+ void MakeNewFaces( TElemIntPairVec& theNew2OldFaces,
+ TNodeIntPairVec& theNew2OldNodes,
+ const double theSign,
+ const bool theOptimize );
+
+ void IntersectNewEdges( const CutFace& theCFace );
private:
bool & isCollinear );
bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint );
bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes );
- void intersectNewEdges( const CutFace& theCFace );
const SMDS_MeshNode* createNode( const gp_XYZ& p );
};
*/
//================================================================================
- const SMDS_MeshNode* Intersector::createNode( const gp_XYZ& p )
+ const SMDS_MeshNode* Intersector::Algo::createNode( const gp_XYZ& p )
{
const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() );
n->setIsMarked( true ); // cut nodes are marked
*/
//================================================================================
- void Intersector::addLink( CutLink& link )
+ void Intersector::Algo::addLink( CutLink& link )
{
link.myIndex = 0;
const CutLink* added = & myCutLinks.Added( link );
*/
//================================================================================
- bool Intersector::findLink( CutLink& link )
+ bool Intersector::Algo::findLink( CutLink& link )
{
link.myIndex = 0;
while ( myCutLinks.Contains( link ))
*/
//================================================================================
- bool Intersector::isPlaneIntersected( const gp_XYZ& n2,
- const double d2,
- const std::vector< SMESH_NodeXYZ >& nodes1,
- std::vector< double > & dist1,
- int & nbOnPlane1,
- int & iNotOnPlane1)
+ bool Intersector::Algo::isPlaneIntersected( const gp_XYZ& n2,
+ const double d2,
+ const std::vector< SMESH_NodeXYZ >& nodes1,
+ std::vector< double > & dist1,
+ int & nbOnPlane1,
+ int & iNotOnPlane1)
{
iNotOnPlane1 = nbOnPlane1 = 0;
dist1.resize( nodes1.size() );
*/
//================================================================================
- void Intersector::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
- const std::vector< double >& dist,
- const int nbOnPln,
- const int iMaxCoo,
- double * u,
- int* iE)
+ void Intersector::Algo::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes,
+ const std::vector< double >& dist,
+ const int nbOnPln,
+ const int iMaxCoo,
+ double * u,
+ int* iE)
{
if ( nbOnPln == 3 )
{
*/
//================================================================================
- void Intersector::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
- const std::vector< double > & dist,
- CutLink& link )
+ void Intersector::Algo::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes,
+ const std::vector< double > & dist,
+ CutLink& link )
{
int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 );
CutLink link2 = link;
*/
//================================================================================
- void Intersector::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
- const std::vector< double > & dist1,
- const int iEdge1,
- const SMDS_MeshElement* face2,
- CutLink& link1)
+ void Intersector::Algo::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1,
+ const std::vector< double > & dist1,
+ const int iEdge1,
+ const SMDS_MeshElement* face2,
+ CutLink& link1)
{
const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size();
const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ];
*/
//================================================================================
- void Intersector::replaceIntNode( const SMDS_MeshNode* nToKeep,
- const SMDS_MeshNode* nToRemove )
+ void Intersector::Algo::replaceIntNode( const SMDS_MeshNode* nToKeep,
+ const SMDS_MeshNode* nToRemove )
{
if ( nToKeep == nToRemove )
return;
if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID
- myRemove2KeepNodes.Bind((void*) nToKeep, nToRemove );
+ myRemove2KeepNodes.Bind( nToKeep, nToRemove );
else
- myRemove2KeepNodes.Bind((void*) nToRemove, nToKeep );
+ myRemove2KeepNodes.Bind( nToRemove, nToKeep );
}
//================================================================================
*/
//================================================================================
- void Intersector::computeIntPoint( const double u1,
- const double u2,
- const int iE1,
- const int iE2,
- CutLink & link,
- const SMDS_MeshNode* & node1,
- const SMDS_MeshNode* & node2)
+ void Intersector::Algo::computeIntPoint( const double u1,
+ const double u2,
+ const int iE1,
+ const int iE2,
+ CutLink & link,
+ const SMDS_MeshNode* & node1,
+ const SMDS_MeshNode* & node2)
{
if ( u1 > u2 + myTol )
{
*/
//================================================================================
- void Intersector::cutCollinearLink( const int iNotOnPlane1,
- const std::vector< SMESH_NodeXYZ >& nodes1,
- const SMDS_MeshElement* face2,
- const CutLink& link1,
- const CutLink& link2)
+ void Intersector::Algo::cutCollinearLink( const int iNotOnPlane1,
+ const std::vector< SMESH_NodeXYZ >& nodes1,
+ const SMDS_MeshElement* face2,
+ const CutLink& link1,
+ const CutLink& link2)
{
int iN1 = ( iNotOnPlane1 + 1 ) % 3;
*/
//================================================================================
- void Intersector::setPlaneIndices( const gp_XYZ& planeNorm )
+ void Intersector::Algo::setPlaneIndices( const gp_XYZ& planeNorm )
{
switch ( MaxIndex( planeNorm )) {
case 1: myInd1 = 2; myInd2 = 3; break;
*/
//================================================================================
- void Intersector::Cut( const SMDS_MeshElement* face1,
- const SMDS_MeshElement* face2,
- const int nbCommonNodes)
+ void Intersector::Algo::Cut( const SMDS_MeshElement* face1,
+ const SMDS_MeshElement* face2,
+ const int nbCommonNodes)
{
myFace1 = face1;
myFace2 = face2;
if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode;
cf1.AddPoint( link1, link2, myTol );
+ if ( l1n1 ) link1.Set( l1n1, l1n2, face2 );
+ if ( l2n1 ) link2.Set( l2n1, l2n2, face2 );
cf2.AddPoint( link1, link2, myTol );
}
else
return;
}
+ //================================================================================
+ /*!
+ * \brief Store a face cut by a line given by its ends
+ * accompanied by indices of intersected face edges.
+ * Edge index is <0 if a line end is inside the face.
+ * \param [in] face - a face to cut
+ * \param [inout] lineEnd1 - line end coordinates + optional node existing at this point
+ * \param [in] edgeIndex1 - index of face edge cut by lineEnd1
+ * \param [inout] lineEnd2 - line end coordinates + optional node existing at this point
+ * \param [in] edgeIndex2 - index of face edge cut by lineEnd2
+ */
+ //================================================================================
+
+ void Intersector::Algo::Cut( const SMDS_MeshElement* face,
+ SMESH_NodeXYZ& lineEnd1,
+ int edgeIndex1,
+ SMESH_NodeXYZ& lineEnd2,
+ int edgeIndex2 )
+ {
+ if ( lineEnd1.Node() && lineEnd2.Node() &&
+ face->GetNodeIndex( lineEnd1.Node() ) >= 0 &&
+ face->GetNodeIndex( lineEnd2.Node() ) >= 0 )
+ return; // intersection at a face node or edge
+
+ if ((int) myNormals.size() <= face->GetID() )
+ const_cast< std::vector< gp_XYZ >& >( myNormals ).resize( face->GetID() + 1 );
+
+ const CutFace& cf = myCutFaces.Added( CutFace( face ));
+ cf.InitLinks();
+
+ // look for intersection nodes coincident with line ends
+ CutLink links[2];
+ for ( int is2nd = 0; is2nd < 2; ++is2nd )
+ {
+ SMESH_NodeXYZ& lineEnd = is2nd ? lineEnd2 : lineEnd1;
+ int edgeIndex = is2nd ? edgeIndex2 : edgeIndex1;
+ CutLink & link = links[ is2nd ];
+
+ link.myIntNode = lineEnd;
+
+ for ( size_t i = ( edgeIndex < 0 ? 3 : 0 ); i < cf.myLinks.size(); ++i )
+ if ( coincide( lineEnd, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), myTol ))
+ {
+ link.myIntNode = cf.myLinks[i].myNode1;
+ break;
+ }
+
+ if ( edgeIndex >= 0 )
+ {
+ link.Set( face->GetNode ( edgeIndex ),
+ face->GetNodeWrap( edgeIndex + 1 ),
+ /*cuttingFace=*/0);
+ findLink( link );
+ }
+
+ if ( !link.myIntNode )
+ link.myIntNode.Set( createNode( lineEnd ));
+
+ lineEnd._node = link.IntNode();
+
+ if ( link.myNode[0] )
+ addLink( link );
+ }
+
+ cf.AddEdge( links[0], links[1], /*face=*/0, /*nbOnPlane=*/0, /*iNotOnPlane=*/-1 );
+ }
+
//================================================================================
/*!
* \brief Intersect two 2D line segments
*/
//================================================================================
- bool Intersector::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
- const gp_XY s2p0, const gp_XY s2p1,
- double & t1, double & t2,
- bool & isCollinear )
+ bool Intersector::Algo::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1,
+ const gp_XY s2p0, const gp_XY s2p1,
+ double & t1, double & t2,
+ bool & isCollinear )
{
gp_XY u = s1p1 - s1p0;
gp_XY v = s2p1 - s2p0;
*/
//================================================================================
- bool Intersector::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint )
+ bool Intersector::Algo::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint )
{
int i01 = iE1, i11 = ( iE1 + 1 ) % 3;
int i02 = iE2, i12 = ( iE2 + 1 ) % 3;
*/
//================================================================================
- bool Intersector::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes )
+ bool Intersector::Algo::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes )
{
double bc1, bc2;
SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ),
p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ),
bc1, bc2 );
- return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
+ //return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. );
+ return ( myTol < bc1 && myTol < bc2 && bc1 + bc2 + myTol < 1. );
}
//================================================================================
*/
//================================================================================
- void Intersector::cutCoplanar()
+ void Intersector::Algo::cutCoplanar()
{
// find intersections of edges
}
return;
- } // Intersector::cutCoplanar()
+ } // Intersector::Algo::cutCoplanar()
//================================================================================
/*!
*/
//================================================================================
- void Intersector::intersectNewEdges( const CutFace& cf )
+ void Intersector::Algo::IntersectNewEdges( const CutFace& cf )
{
IntPoint2D intPoint;
if ( cf.NbInternalEdges() < 2 )
return;
+ if ( myNodes1.empty() )
+ {
+ myNodes1.resize(2);
+ myNodes2.resize(2);
+ }
+
const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ];
setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane
size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2;
- for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 )
+ size_t i1 = 3;
+ while ( cf.myLinks[i1-1].IsInternal() && i1 > 0 )
+ --i1;
+
+ for ( ; i1 < cf.myLinks.size(); ++i1 )
{
if ( !cf.myLinks[i1].IsInternal() )
continue;
}
}
if ( cf.myLinks.size() >= limit )
- throw SALOME_Exception( "Infinite loop in Intersector::intersectNewEdges()" );
+ throw SALOME_Exception( "Infinite loop in Intersector::Algo::IntersectNewEdges()" );
}
++i1; // each internal edge encounters twice
}
*/
//================================================================================
- void Intersector::MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces,
- SMESH_MeshAlgos::TNPairVec& theNew2OldNodes,
- const double theSign)
+ void Intersector::Algo::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
+ SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
+ const double theSign,
+ const bool theOptimize)
{
+ // fill theNew2OldFaces if empty
+ TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin();
+ if ( theNew2OldFaces.empty() )
+ for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+ {
+ const CutFace& cf = *cutFacesIt;
+ int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
+ if ((int) theNew2OldFaces.size() <= index )
+ theNew2OldFaces.resize( index + 1 );
+ theNew2OldFaces[ index ] = std::make_pair( cf.myInitFace, index );
+ }
+
// unmark all nodes except intersection ones
for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); )
// intersect edges added to myCutFaces
- TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin();
- for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+ for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
{
const CutFace& cf = *cutFacesIt;
cf.ReplaceNodes( myRemove2KeepNodes );
- intersectNewEdges( cf );
+ IntersectNewEdges( cf );
}
// make new faces
EdgeLoopSet loopSet;
- SMESH_MeshAlgos::Triangulate triangulator;
+ SMESH_MeshAlgos::Triangulate triangulator( theOptimize );
std::vector< EdgePart > cutOffLinks;
TLinkMap cutOffCoplanarLinks;
std::vector< const CutFace* > touchedFaces;
- SMESH_MeshAlgos::TEPairVec::value_type new2OldTria;
+ SMESH_MeshAlgos::TElemIntPairVec::value_type new2OldTria;
CutFace cutFace(0);
std::vector< const SMDS_MeshNode* > nodes;
std::vector<const SMDS_MeshElement *> faces;
// avoid loops that are not connected to boundary edges of cf.myInitFace
if ( cf.RemoveInternalLoops( loopSet ))
{
- intersectNewEdges( cf );
+ IntersectNewEdges( cf );
cf.MakeLoops( loopSet, normal );
}
// erase loops that are cut off by face intersections
// remove split faces
for ( size_t id = 1; id < theNew2OldFaces.size(); ++id )
{
- if ( theNew2OldFaces[id].first )
+ if ( theNew2OldFaces[id].first ||
+ theNew2OldFaces[id].second == 0 )
continue;
if ( const SMDS_MeshElement* f = myMesh->FindElement( id ))
myMesh->RemoveFreeElement( f );
}
- // remove face connected to cut off parts of cf.myInitFace
+ // remove faces that are merged off
+ for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt )
+ {
+ const CutFace& cf = *cutFacesIt;
+ if ( !cf.myLinks.empty() || cf.myInitFace->IsNull() )
+ continue;
+
+ nodes.assign( cf.myInitFace->begin_nodes(), cf.myInitFace->end_nodes() );
+ for ( size_t i = 0; i < nodes.size(); ++i )
+ {
+ const SMDS_MeshNode* n = nodes[ i ];
+ while ( myRemove2KeepNodes.IsBound( n ))
+ n = myRemove2KeepNodes( n );
+ if ( n != nodes[ i ] && cf.myInitFace->GetNodeIndex( n ) >= 0 )
+ {
+ theNew2OldFaces[ cf.myInitFace->GetID() ].first = 0;
+ myMesh->RemoveFreeElement( cf.myInitFace );
+ break;
+ }
+ }
+ }
+
+ // remove faces connected to cut off parts of cf.myInitFace
nodes.resize(2);
for ( size_t i = 0; i < cutOffLinks.size(); ++i )
if ( nodes[0] != nodes[1] &&
myMesh->GetElementsByNodes( nodes, faces ))
{
- if ( cutOffLinks[i].myFace &&
- cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
- faces.size() == 2 )
+ if ( // cutOffLinks[i].myFace &&
+ cutOffLinks[i].myIndex != EdgePart::_COPLANAR &&
+ faces.size() != 1 )
continue;
for ( size_t iF = 0; iF < faces.size(); ++iF )
{
for ( size_t i = 0; i < touchedFaces.size(); ++i )
{
const CutFace& cf = *touchedFaces[i];
+ if ( cf.myInitFace->IsNull() )
+ continue;
int index = cf.myInitFace->GetID(); // index in theNew2OldFaces
if ( !theNew2OldFaces[ index ].first )
continue; // already cut off
+ cf.InitLinks();
if ( !cf.ReplaceNodes( myRemove2KeepNodes ))
- continue; // just keep as is
+ {
+ if ( cf.myLinks.size() == 3 &&
+ cf.myInitFace->GetNodeIndex( cf.myLinks[0].myNode1 ) >= 0 &&
+ cf.myInitFace->GetNodeIndex( cf.myLinks[1].myNode1 ) >= 0 &&
+ cf.myInitFace->GetNodeIndex( cf.myLinks[2].myNode1 ) >= 0 )
+ continue; // just keep as is
+ }
if ( cf.myLinks.size() == 3 )
{
// add used new nodes to theNew2OldNodes
- SMESH_MeshAlgos::TNPairVec::value_type new2OldNode;
- new2OldNode.second = NULL;
+ SMESH_MeshAlgos::TNodeIntPairVec::value_type new2OldNode;
+ new2OldNode.second = 0;
for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt )
{
const CutLink& link = *cutLinksIt;
return;
}
+ //================================================================================
+ Intersector::Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals )
+ {
+ myAlgo = new Algo( mesh, tol, normals );
+ }
+ //================================================================================
+ Intersector::~Intersector()
+ {
+ delete myAlgo;
+ }
+ //================================================================================
+ //! compute cut of two faces of the mesh
+ void Intersector::Cut( const SMDS_MeshElement* face1,
+ const SMDS_MeshElement* face2,
+ const int nbCommonNodes )
+ {
+ myAlgo->Cut( face1, face2, nbCommonNodes );
+ }
+ //================================================================================
+ //! store a face cut by a line given by its ends
+ // accompanied by indices of intersected face edges.
+ // Edge index is <0 if a line end is inside the face.
+ void Intersector::Cut( const SMDS_MeshElement* face,
+ SMESH_NodeXYZ& lineEnd1,
+ int edgeIndex1,
+ SMESH_NodeXYZ& lineEnd2,
+ int edgeIndex2 )
+ {
+ myAlgo->Cut( face, lineEnd1, edgeIndex1, lineEnd2, edgeIndex2 );
+ }
+ //================================================================================
+ //! split all face intersected by Cut() methods
+ void Intersector::MakeNewFaces( SMESH_MeshAlgos::TElemIntPairVec& theNew2OldFaces,
+ SMESH_MeshAlgos::TNodeIntPairVec& theNew2OldNodes,
+ const double theSign,
+ const bool theOptimize )
+ {
+ myAlgo->MakeNewFaces( theNew2OldFaces, theNew2OldNodes, theSign, theOptimize );
+ }
+ //================================================================================
+ //! Cut a face by planes, whose normals point to parts to keep
+ bool Intersector::CutByPlanes(const SMDS_MeshElement* theFace,
+ const std::vector< gp_Ax1 > & thePlanes,
+ const double theTol,
+ std::vector< TFace > & theNewFaceConnectivity )
+ {
+ theNewFaceConnectivity.clear();
+
+ // check if theFace is wholly cut off
+ std::vector< SMESH_NodeXYZ > facePoints( theFace->begin_nodes(), theFace->end_nodes() );
+ facePoints.resize( theFace->NbCornerNodes() );
+ for ( size_t iP = 0; iP < thePlanes.size(); ++iP )
+ {
+ size_t nbOut = 0;
+ const gp_Pnt& O = thePlanes[iP].Location();
+ for ( size_t i = 0; i < facePoints.size(); ++i )
+ {
+ gp_Vec Op( O, facePoints[i] );
+ nbOut += ( Op * thePlanes[iP].Direction() <= 0 );
+ }
+ if ( nbOut == facePoints.size() )
+ return true;
+ }
+
+ // copy theFace into a temporary mesh
+ SMDS_Mesh mesh;
+ Bnd_B3d faceBox;
+ std::vector< const SMDS_MeshNode* > faceNodes;
+ faceNodes.resize( facePoints.size() );
+ for ( size_t i = 0; i < facePoints.size(); ++i )
+ {
+ const SMESH_NodeXYZ& n = facePoints[i];
+ faceNodes[i] = mesh.AddNode( n.X(), n.Y(), n.Z() );
+ faceBox.Add( n );
+ }
+ const SMDS_MeshElement* faceToCut = 0;
+ switch ( theFace->NbCornerNodes() )
+ {
+ case 3:
+ faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2] );
+ break;
+ case 4:
+ faceToCut = mesh.AddFace( faceNodes[0], faceNodes[1], faceNodes[2], faceNodes[3] );
+ break;
+ default:
+ faceToCut = mesh.AddPolygonalFace( faceNodes );
+ }
+
+ std::vector< gp_XYZ > normals( 2 + thePlanes.size() );
+ SMESH_MeshAlgos::FaceNormal( faceToCut, normals[ faceToCut->GetID() ]);
+
+ // add faces corresponding to thePlanes
+ std::vector< const SMDS_MeshElement* > planeFaces;
+ double faceSize = Sqrt( faceBox.SquareExtent() );
+ gp_XYZ center = 0.5 * ( faceBox.CornerMin() + faceBox.CornerMax() );
+ for ( size_t i = 0; i < thePlanes.size(); ++i )
+ {
+ gp_Ax2 plnAx( thePlanes[i].Location(), thePlanes[i].Direction() );
+ gp_XYZ O = plnAx.Location().XYZ();
+ gp_XYZ X = plnAx.XDirection().XYZ();
+ gp_XYZ Y = plnAx.YDirection().XYZ();
+ gp_XYZ Z = plnAx.Direction().XYZ();
+
+ double dot = ( O - center ) * Z;
+ gp_XYZ o = center + Z * dot; // center projected to a plane
+
+ gp_XYZ p1 = o + X * faceSize * 2;
+ gp_XYZ p2 = o + Y * faceSize * 2;
+ gp_XYZ p3 = o - (X + Y ) * faceSize * 2;
+
+ const SMDS_MeshNode* n1 = mesh.AddNode( p1.X(), p1.Y(), p1.Z() );
+ const SMDS_MeshNode* n2 = mesh.AddNode( p2.X(), p2.Y(), p2.Z() );
+ const SMDS_MeshNode* n3 = mesh.AddNode( p3.X(), p3.Y(), p3.Z() );
+ planeFaces.push_back( mesh.AddFace( n1, n2, n3 ));
+
+ normals[ planeFaces.back()->GetID() ] = thePlanes[i].Direction().XYZ();
+ }
+
+ // cut theFace
+ Algo algo ( &mesh, theTol, normals );
+ for ( size_t i = 0; i < planeFaces.size(); ++i )
+ {
+ algo.Cut( faceToCut, planeFaces[i], 0 );
+ }
+
+ // retrieve a result
+ SMESH_MeshAlgos::TElemIntPairVec new2OldFaces;
+ SMESH_MeshAlgos::TNodeIntPairVec new2OldNodes;
+ TCutFaceMap::const_iterator cutFacesIt= algo.myCutFaces.cbegin();
+ for ( ; cutFacesIt != algo.myCutFaces.cend(); ++cutFacesIt )
+ {
+ const CutFace& cf = *cutFacesIt;
+ if ( cf.myInitFace != faceToCut )
+ continue;
+
+ if ( !cf.IsCut() )
+ {
+ theNewFaceConnectivity.push_back( facePoints );
+ break;
+ }
+
+ // intersect cut lines
+ algo.IntersectNewEdges( cf );
+
+ // form loops of new faces
+ EdgeLoopSet loopSet;
+ cf.MakeLoops( loopSet, normals[ faceToCut->GetID() ]);
+
+ // erase loops that are cut off by thePlanes
+ const double sign = 1;
+ std::vector< EdgePart > cutOffLinks;
+ TLinkMap cutOffCoplanarLinks;
+ cf.CutOffLoops( loopSet, sign, normals, cutOffLinks, cutOffCoplanarLinks );
+
+ for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL )
+ {
+ EdgeLoop& loop = loopSet.myLoops[ iL ];
+ if ( loop.myLinks.size() > 0 )
+ {
+ facePoints.clear();
+ for ( SMDS_NodeIteratorPtr nIt = loop.nodeIterator(); nIt->more(); )
+ {
+ const SMDS_MeshNode* n = nIt->next();
+ facePoints.push_back( n );
+ int iN = faceToCut->GetNodeIndex( n );
+ if ( iN < 0 )
+ facePoints.back()._node = 0; // an intersection point
+ else
+ facePoints.back()._node = theFace->GetNode( iN );
+ }
+ theNewFaceConnectivity.push_back( facePoints );
+ }
+ }
+ break;
+ }
+
+ return theNewFaceConnectivity.empty();
+ }
+
+} // namespace SMESH_MeshAlgos
+
+namespace
+{
//================================================================================
/*!
* \brief Debug
// if ( !myLinks[i].IsInternal() )
// myLinks[ i ].myFace = cutterFace;
// else
- myLinks[ i ].ReplaceCoplanar( newEdge );
- myLinks[ i+1 ].ReplaceCoplanar( newEdge );
+ myLinks[ i ].ReplaceCoplanar( newEdge );
+ if ( myLinks[i].IsInternal() && i+1 < myLinks.size() )
+ myLinks[ i+1 ].ReplaceCoplanar( newEdge );
return;
}
i += myLinks[i].IsInternal();
bool replaced = false;
for ( size_t i = 0; i < myLinks.size(); ++i )
{
- while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode1 ))
- replaced = ( myLinks[i].myNode1 = theRm2KeepMap((Standard_Address) myLinks[i].myNode1 ));
+ while ( theRm2KeepMap.IsBound( myLinks[i].myNode1 ))
+ replaced = ( myLinks[i].myNode1 = theRm2KeepMap( myLinks[i].myNode1 ));
- while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode2 ))
- replaced = ( myLinks[i].myNode2 = theRm2KeepMap((Standard_Address) myLinks[i].myNode2 ));
+ while ( theRm2KeepMap.IsBound( myLinks[i].myNode2 ))
+ replaced = ( myLinks[i].myNode2 = theRm2KeepMap( myLinks[i].myNode2 ));
}
//if ( replaced ) // remove equal links
//================================================================================
/*!
* \brief Remove loops that are not connected to boundary edges of myFace by
- * adding edges connecting these loops to the boundary
+ * adding edges connecting these loops to the boundary.
+ * Such loops must be removed as they form polygons with ring topology.
*/
//================================================================================
while ( prevNbReached < nbReachedLoops );
- // add links connecting internal loops with the boundary ones
for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
{
EdgeLoop& loop = theLoops.myLoops[ iL ];
- if ( loop.myIsBndConnected )
+ if ( loop.myIsBndConnected || loop.myLinks.size() == 0 )
continue;
+ if ( loop.myHasPending )
+ {
+ // try to join the loop to another one, with which it contacts at a node
+
+ // look for a node where the loop reverses
+ const EdgePart* edgePrev = loop.myLinks.back();
+ for ( size_t iE = 0; iE < loop.myLinks.size(); edgePrev = loop.myLinks[ iE++ ] )
+ {
+ if ( !edgePrev->IsTwin( *loop.myLinks[ iE ]))
+ continue;
+ const SMDS_MeshNode* reverseNode = edgePrev->myNode2;
+
+ // look for a loop including reverseNode
+ size_t iContactEdge2; // index(+1) of edge starting at reverseNode
+ for ( size_t iL2 = 0; iL2 < theLoops.myNbLoops; ++iL2 )
+ {
+ if ( iL == iL2 )
+ continue;
+ EdgeLoop& loop2 = theLoops.myLoops[ iL2 ];
+ if ( ! ( iContactEdge2 = loop2.Contains( reverseNode )))
+ continue;
+
+ // insert loop2 into the loop
+ theLoops.Join( loop, iE, loop2, iContactEdge2 - 1 );
+ break;
+ }
+ if ( loop.myIsBndConnected )
+ break;
+ }
+
+ if ( loop.myIsBndConnected )
+ continue;
+ }
+
+ // add links connecting internal loops with the boundary ones
+
// find a pair of closest nodes
const SMDS_MeshNode *closestNode1, *closestNode2;
double minDist = 1e100;
{
theLoops.AddNewLoop();
theLoops.AddEdge( myLinks[0] );
- theLoops.AddEdge( myLinks[1] );
- theLoops.AddEdge( myLinks[2] );
+ if ( myLinks[0].myNode2 == myLinks[1].myNode1 )
+ {
+ theLoops.AddEdge( myLinks[1] );
+ theLoops.AddEdge( myLinks[2] );
+ }
+ else
+ {
+ theLoops.AddEdge( myLinks[2] );
+ theLoops.AddEdge( myLinks[1] );
+ }
return;
}
while ( !theLoops.AllEdgesUsed() )
{
- theLoops.AddNewLoop();
+ EdgeLoop& loop = theLoops.AddNewLoop();
// add 1st edge to a new loop
size_t i1;
// choose among candidates
if ( theLoops.myCandidates.size() == 0 )
{
- theLoops.GetLoopOf( lastEdge )->myHasPending = true;
+ loop.myHasPending = bool( twinEdge );
lastEdge = twinEdge;
}
else if ( theLoops.myCandidates.size() == 1 )
}
while ( lastNode != firstNode );
+
+ if ( twinEdge == & myLinks[ i1 ])
+ loop.myHasPending = true;
+
} // while ( !theLoops.AllEdgesUsed() )
return;
TLinkMap& theCutOffCoplanarLinks) const
{
EdgePart sideEdge;
+ boost::container::flat_set< const SMDS_MeshElement* > checkedCoplanar;
+
for ( size_t i = 0; i < myLinks.size(); ++i )
{
if ( !myLinks[i].myFace )
loop = theLoops.GetLoopOf( twin );
toErase = ( loop && !loop->myLinks.empty() );
}
+
+ if ( toErase ) // do not erase if cutFace is connected to a co-planar cutFace
+ {
+ checkedCoplanar.clear();
+ for ( size_t iE = 0; iE < myLinks.size() && toErase; ++iE )
+ {
+ if ( !myLinks[iE].myFace || myLinks[iE].myIndex != EdgePart::_COPLANAR )
+ continue;
+ bool isAdded = checkedCoplanar.insert( myLinks[iE].myFace ).second;
+ if ( !isAdded )
+ continue;
+ toErase = ( SMESH_MeshAlgos::NbCommonNodes( myLinks[i ].myFace,
+ myLinks[iE].myFace ) < 1 );
+ }
+ }
}
if ( toErase )
{
if ( !loop->myLinks[ iE ]->myFace &&
!loop->myLinks[ iE ]->IsInternal() )// &&
- // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
- // !loop->myLinks[ iE ]->myNode2->isMarked() )
+ // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked
+ // !loop->myLinks[ iE ]->myNode2->isMarked() )
{
int i = loop->myLinks[ iE ]->myIndex;
sideEdge.Set( myInitFace->GetNode ( i ),
//================================================================================
/*!
- * \brief Replace _COPLANAR cut edge by _INTERNAL oe vice versa
+ * \brief Replace _COPLANAR cut edge by _INTERNAL or vice versa
*/
//================================================================================
if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL )
{
//check if the faces are connected
- int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size();
+ int nbCommonNodes = 0;
+ if ( e.myFace && myFace )
+ nbCommonNodes = SMESH_MeshAlgos::NbCommonNodes( e.myFace, myFace );
bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) ||
( myIndex == _COPLANAR && nbCommonNodes < 2 ));
if ( toReplace )
/*!
* \brief Create an offsetMesh of given faces
* \param [in] faceIt - the input faces
- * \param [out] new2OldFaces - history of faces
- * \param [out] new2OldNodes - history of nodes
+ * \param [out] new2OldFaces - history of faces (new face -> old face ID)
+ * \param [out] new2OldNodes - history of nodes (new node -> old node ID)
* \return SMDS_Mesh* - the new offset mesh, a caller should delete
*/
//================================================================================
SMDS_Mesh& theSrcMesh,
const double theOffset,
const bool theFixIntersections,
- TEPairVec& theNew2OldFaces,
- TNPairVec& theNew2OldNodes)
+ TElemIntPairVec& theNew2OldFaces,
+ TNodeIntPairVec& theNew2OldNodes)
{
- SMDS_Mesh* newMesh = new SMDS_Mesh;
- theNew2OldFaces.clear();
- theNew2OldNodes.clear();
- theNew2OldFaces.push_back
- ( std::make_pair(( const SMDS_MeshElement*) 0,
- ( const SMDS_MeshElement*) 0)); // to have index == face->GetID()
-
if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 )
throw SALOME_Exception( "Offset of quadratic mesh not supported" );
if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() )
throw SALOME_Exception( "Offset of non-triangular mesh not supported" );
+ SMDS_Mesh* newMesh = new SMDS_Mesh;
+ theNew2OldFaces.clear();
+ theNew2OldNodes.clear();
+ theNew2OldFaces.push_back
+ ( std::make_pair(( const SMDS_MeshElement*) 0, 0)); // to have index == face->GetID()
+
// copy input faces to the newMesh keeping IDs of nodes
double minNodeDist = 1e100;
{
SMESH_NodeXYZ xyz( nodes[i] );
newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() );
- theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i] ));
+ theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i]->GetID() ));
nodes[i] = newNode;
}
}
default:
continue;
}
- theNew2OldFaces.push_back( std::make_pair( newFace, face ));
+ theNew2OldFaces.push_back( std::make_pair( newFace, face->GetID() ));
SMESH_NodeXYZ pPrev = nodes.back(), p;
for ( size_t i = 0; i < nodes.size(); ++i )
{
p.Set( nodes[i] );
double dist = ( pPrev - p ).SquareModulus();
- if ( dist > std::numeric_limits<double>::min() )
+ if ( dist < minNodeDist && dist > std::numeric_limits<double>::min() )
minNodeDist = dist;
pPrev = p;
}
std::vector< gp_XYZ > normals( theNew2OldFaces.size() );
for ( size_t i = 1; i < normals.size(); ++i )
{
- if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].second, normals[i] ))
+ if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].first, normals[i] ))
normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors
}
- const double tol = 1e-3 * Sqrt( minNodeDist );
const double sign = ( theOffset < 0 ? -1 : +1 );
+ const double tol = Min( 1e-3 * Sqrt( minNodeDist ),
+ 1e-2 * theOffset * sign );
// translate new nodes by normal to input faces
gp_XYZ newXYZ;
{
newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() );
newNode->setIsMarked( true );
- theNew2OldNodes.push_back( std::make_pair( newNode, theNew2OldNodes[i].second ));
+ theNew2OldNodes.push_back( std::make_pair( newNode, 0 ));
multiPos.emplace_back( newNode );
}
}
// mark all new nodes located closer than theOffset from theSrcMesh
}
+ removeSmallFaces( newMesh, theNew2OldFaces, tol*tol );
+
// ==================================================
// find self-intersections of new faces and fix them
// ==================================================
Intersector intersector( newMesh, tol, normals );
std::vector< const SMDS_MeshElement* > closeFaces;
- std::vector< const SMDS_MeshNode* > faceNodes;
+ std::vector< SMESH_NodeXYZ > faceNodes;
Bnd_B3d faceBox;
+
for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF )
{
const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first;
closeFaces.clear();
faceBox.Clear();
for ( size_t i = 0; i < faceNodes.size(); ++i )
- faceBox.Add( SMESH_NodeXYZ( faceNodes[i] ));
+ faceBox.Add( faceNodes[i] );
faceBox.Enlarge( tol );
fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces );
// do not intersect connected faces if they have no concave nodes
int nbCommonNodes = 0;
for ( size_t iN = 0; iN < faceNodes.size(); ++iN )
- nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN] ) >= 0 );
+ nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN].Node() ) >= 0 );
if ( !isConcaveNode1 )
{
break;
if ( !isConcaveNode2 && nbCommonNodes > 0 )
- continue;
+ {
+ if ( normals[ newFace->GetID() ] * normals[ closeFace->GetID() ] < 1.0 )
+ continue; // not co-planar
+ }
}
intersector.Cut( newFace, closeFace, nbCommonNodes );
}
}
- intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign );
+ intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign, /*optimize=*/true );
return newMesh;
}
