-// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012 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
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
// SMESH SMESH : idl implementation based on 'SMESH' unit's classes
// File : SMESH_MeshEditor.cxx
#include <algorithm>
#include <sstream>
+#include <Standard_Failure.hxx>
+#include <Standard_ErrorHandler.hxx>
+
#define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
using namespace std;
return false;
}
+//================================================================================
+/*!
+ * \brief Reorient faces.
+ * \param theFaces - the faces to reorient. If empty the whole mesh is meant
+ * \param theDirection - desired direction of normal of \a theFace
+ * \param theFace - one of \a theFaces that sould be orientated according to
+ * \a theDirection and whose orientation defines orientation of other faces
+ * \return number of reoriented faces.
+ */
+//================================================================================
+
+int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces,
+ const gp_Dir& theDirection,
+ const SMDS_MeshElement * theFace)
+{
+ int nbReori = 0;
+ if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori;
+
+ if ( theFaces.empty() )
+ {
+ SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true);
+ while ( fIt->more() )
+ theFaces.insert( theFaces.end(), fIt->next() );
+ }
+
+ // orient theFace according to theDirection
+ gp_XYZ normal;
+ SMESH_Algo::FaceNormal( theFace, normal, /*normalized=*/false );
+ if ( normal * theDirection.XYZ() < 0 )
+ nbReori += Reorient( theFace );
+
+ // Orient other faces
+
+ set< const SMDS_MeshElement* > startFaces;
+ TIDSortedElemSet avoidSet;
+ set< SMESH_TLink > checkedLinks;
+ pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew;
+
+ if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces
+ theFaces.erase( theFace );
+ startFaces.insert( theFace );
+
+ set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin();
+ while ( startFace != startFaces.end() )
+ {
+ theFace = *startFace;
+ const int nbNodes = theFace->NbCornerNodes();
+
+ avoidSet.clear();
+ avoidSet.insert(theFace);
+
+ NLink link( theFace->GetNode( 0 ), 0 );
+ for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace
+ {
+ link.second = theFace->GetNode(( i+1 ) % nbNodes );
+ linkIt_isNew = checkedLinks.insert( link );
+ if ( !linkIt_isNew.second )
+ {
+ // link has already been checked and won't be encountered more
+ // if the group (theFaces) is manifold
+ checkedLinks.erase( linkIt_isNew.first );
+ }
+ else
+ {
+ int nodeInd1, nodeInd2;
+ const SMDS_MeshElement* otherFace = FindFaceInSet( link.first, link.second,
+ theFaces, avoidSet,
+ & nodeInd1, & nodeInd2);
+ if ( otherFace && otherFace != theFace)
+ {
+ // link must be reversed in otherFace if orientation ot otherFace
+ // is same as that of theFace
+ if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
+ {
+ // cout << "Reorient " << otherFace->GetID() << " near theFace=" <<theFace->GetID()
+ // << " \tlink( " << link.first->GetID() << " " << link.second->GetID() << endl;
+ nbReori += Reorient( otherFace );
+ }
+ startFaces.insert( otherFace );
+ if ( theFaces.size() > 1 ) // leave 1 face to prevent finding not selected faces
+ theFaces.erase( otherFace );
+ }
+ }
+ std::swap( link.first, link.second );
+ }
+ startFaces.erase( startFace );
+ startFace = startFaces.begin();
+ }
+ return nbReori;
+}
+
//=======================================================================
//function : getBadRate
//purpose :
break;
}
case SMDSEntity_Quad_Quadrangle: { // sweep Quadratic QUADRANGLE --->
- if ( nbDouble != 4 ) break;
if( nbSame == 0 ) {
+ if ( nbDouble != 4 ) break;
// ---> hexahedron with 20 nodes
aNewElem = aMesh->AddVolume (prevNod[0], prevNod[1], prevNod[2], prevNod[3],
nextNod[0], nextNod[1], nextNod[2], nextNod[3],
return;
}
else if( nbSame == 2 ) {
+ if ( nbDouble != 2 ) break;
// ---> 2d order Pentahedron with 15 nodes
int n1,n2,n4,n5;
if ( prevNod[ iBeforeSame ] == nextNod[ iBeforeSame ] ) {
const SMDS_MeshElement* currentElem = NULL;
int totalNbEdges = theTrack->NbEdges();
SMDS_ElemIteratorPtr nIt;
- bool isClosed = false;
//check start node
if( !theTrack->GetMeshDS()->Contains(theN1) ) {
//case of the closed mesh
if(currentNode == theN1) {
nbEdges++;
- isClosed = true;
break;
}
{
public:
- ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius );
- void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems);
+ ElementBndBoxTree(const SMDS_Mesh& mesh,
+ SMDSAbs_ElementType elemType,
+ SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(),
+ double tolerance = NodeRadius );
+ void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems );
void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems);
+ void getElementsInSphere ( const gp_XYZ& center,
+ const double radius, TIDSortedElemSet& foundElems);
+ size_t getSize() { return std::max( _size, _elements.size() ); }
~ElementBndBoxTree();
protected:
- ElementBndBoxTree() {}
+ ElementBndBoxTree():_size(0) {}
SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; }
- void buildChildrenData();
- Bnd_B3d* buildRootBox();
+ void buildChildrenData();
+ Bnd_B3d* buildRootBox();
private:
//!< Bounding box of element
struct ElementBox : public Bnd_B3d
ElementBox(const SMDS_MeshElement* elem, double tolerance);
};
vector< ElementBox* > _elements;
+ size_t _size;
};
//================================================================================
while ( elemIt->more() )
_elements.push_back( new ElementBox( elemIt->next(),tolerance ));
- if ( _elements.size() > MaxNbElemsInLeaf )
- compute();
- else
- myIsLeaf = true;
+ compute();
}
//================================================================================
}
_elements[i]->_refCount--;
}
- _elements.clear();
+ _size = _elements.size();
+ SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory
for (int j = 0; j < 8; j++)
{
child->myIsLeaf = true;
if ( child->_elements.capacity() - child->_elements.size() > 1000 )
- child->_elements.resize( child->_elements.size() ); // compact
+ SMESHUtils::CompactVector( child->_elements );
}
}
void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point,
TIDSortedElemSet& foundElems)
{
- if ( level() && getBox().IsOut( point.XYZ() ))
+ if ( getBox().IsOut( point.XYZ() ))
return;
if ( isLeaf() )
void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line,
TIDSortedElemSet& foundElems)
{
- if ( level() && getBox().IsOut( line ))
+ if ( getBox().IsOut( line ))
return;
if ( isLeaf() )
}
}
+ //================================================================================
+ /*!
+ * \brief Return elements from leaves intersecting the sphere
+ */
+ //================================================================================
+
+ void ElementBndBoxTree::getElementsInSphere ( const gp_XYZ& center,
+ const double radius,
+ TIDSortedElemSet& foundElems)
+ {
+ if ( getBox().IsOut( center, radius ))
+ return;
+
+ if ( isLeaf() )
+ {
+ for ( int i = 0; i < _elements.size(); ++i )
+ if ( !_elements[i]->IsOut( center, radius ))
+ foundElems.insert( _elements[i]->_element );
+ }
+ else
+ {
+ for (int i = 0; i < 8; i++)
+ ((ElementBndBoxTree*) myChildren[i])->getElementsInSphere( center, radius, foundElems );
+ }
+ }
+
//================================================================================
/*!
* \brief Construct the element box
_refCount = 1;
SMDS_ElemIteratorPtr nIt = elem->nodesIterator();
while ( nIt->more() )
- Add( SMESH_TNodeXYZ( cast2Node( nIt->next() )));
+ Add( SMESH_TNodeXYZ( nIt->next() ));
Enlarge( tolerance );
}
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElements);
virtual TopAbs_State GetPointState(const gp_Pnt& point);
+ virtual const SMDS_MeshElement* FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type );
void GetElementsNearLine( const gp_Ax1& line,
SMDSAbs_ElementType type,
_ebbTree->getElementsNearPoint( point, suspectElems );
TIDSortedElemSet::iterator elem = suspectElems.begin();
for ( ; elem != suspectElems.end(); ++elem )
- if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance ))
+ if ( !SMESH_MeshEditor::IsOut( *elem, point, tolerance ))
foundElements.push_back( *elem );
}
return foundElements.size();
}
+//=======================================================================
+/*!
+ * \brief Find an element of given type most close to the given point
+ *
+ * WARNING: Only face search is implemeneted so far
+ */
+//=======================================================================
+
+const SMDS_MeshElement*
+SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point,
+ SMDSAbs_ElementType type )
+{
+ const SMDS_MeshElement* closestElem = 0;
+
+ if ( type == SMDSAbs_Face )
+ {
+ if ( !_ebbTree || _elementType != type )
+ {
+ if ( _ebbTree ) delete _ebbTree;
+ _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt );
+ }
+ TIDSortedElemSet suspectElems;
+ _ebbTree->getElementsNearPoint( point, suspectElems );
+
+ if ( suspectElems.empty() && _ebbTree->maxSize() > 0 )
+ {
+ gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox().CornerMin() +
+ _ebbTree->getBox().CornerMax() );
+ double radius;
+ if ( _ebbTree->getBox().IsOut( point.XYZ() ))
+ radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize();
+ else
+ radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2;
+ while ( suspectElems.empty() )
+ {
+ _ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems );
+ radius *= 1.1;
+ }
+ }
+ double minDist = std::numeric_limits<double>::max();
+ multimap< double, const SMDS_MeshElement* > dist2face;
+ TIDSortedElemSet::iterator elem = suspectElems.begin();
+ for ( ; elem != suspectElems.end(); ++elem )
+ {
+ double dist = SMESH_MeshEditor::GetDistance( dynamic_cast<const SMDS_MeshFace*>(*elem),
+ point );
+ if ( dist < minDist + 1e-10)
+ {
+ minDist = dist;
+ dist2face.insert( dist2face.begin(), make_pair( dist, *elem ));
+ }
+ }
+ if ( !dist2face.empty() )
+ {
+ multimap< double, const SMDS_MeshElement* >::iterator d2f = dist2face.begin();
+ closestElem = d2f->second;
+ // if there are several elements at the same distance, select one
+ // with GC closest to the point
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ double minDistToGC = 0;
+ for ( ++d2f; d2f != dist2face.end() && fabs( d2f->first - minDist ) < 1e-10; ++d2f )
+ {
+ if ( minDistToGC == 0 )
+ {
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator(closestElem->nodesIterator()),
+ TXyzIterator(), gc ) / closestElem->NbNodes();
+ minDistToGC = point.SquareDistance( gc );
+ }
+ gp_XYZ gc(0,0,0);
+ gc = accumulate( TXyzIterator( d2f->second->nodesIterator()),
+ TXyzIterator(), gc ) / d2f->second->NbNodes();
+ double d = point.SquareDistance( gc );
+ if ( d < minDistToGC )
+ {
+ minDistToGC = d;
+ closestElem = d2f->second;
+ }
+ }
+ // cout << "FindClosestTo( " <<point.X()<<", "<<point.Y()<<", "<<point.Z()<<" ) FACE "
+ // <<closestElem->GetID() << " DIST " << minDist << endl;
+ }
+ }
+ else
+ {
+ // NOT IMPLEMENTED SO FAR
+ }
+ return closestElem;
+}
+
+
//================================================================================
/*!
* \brief Classify the given point in the closed 2D mesh
else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 )
{
gp_Pnt intersectionPoint = intersection.Point(1);
- if ( !SMESH_MeshEditor::isOut( *face, intersectionPoint, tolerance ))
+ if ( !SMESH_MeshEditor::IsOut( *face, intersectionPoint, tolerance ))
u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm )));
}
}
*/
//=======================================================================
-bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
+bool SMESH_MeshEditor::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol )
{
if ( element->GetType() == SMDSAbs_Volume)
{
for ( i = 0; i < nbNodes; ++i )
{
SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] );
- if ( !isOut( &edge, point, tol ))
+ if ( !IsOut( &edge, point, tol ))
return false;
}
return true;
return true;
}
+//=======================================================================
+
+namespace
+{
+ // Position of a point relative to a segment
+ // . .
+ // . LEFT .
+ // . .
+ // VERTEX 1 o----ON-----> VERTEX 2
+ // . .
+ // . RIGHT .
+ // . .
+ enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8,
+ POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX };
+ struct PointPos
+ {
+ PositionName _name;
+ int _index; // index of vertex or segment
+
+ PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {}
+ bool operator < (const PointPos& other ) const
+ {
+ if ( _name == other._name )
+ return ( _index < 0 || other._index < 0 ) ? false : _index < other._index;
+ return _name < other._name;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return of a point relative to a segment
+ * \param point2D - the point to analyze position of
+ * \param xyVec - end points of segments
+ * \param index0 - 0-based index of the first point of segment
+ * \param posToFindOut - flags of positions to detect
+ * \retval PointPos - point position
+ */
+ //================================================================================
+
+ PointPos getPointPosition( const gp_XY& point2D,
+ const gp_XY* segEnds,
+ const int index0 = 0,
+ const int posToFindOut = POS_ALL)
+ {
+ const gp_XY& p1 = segEnds[ index0 ];
+ const gp_XY& p2 = segEnds[ index0+1 ];
+ const gp_XY grad = p2 - p1;
+
+ if ( posToFindOut & POS_VERTEX )
+ {
+ // check if the point2D is at "vertex 1" zone
+ gp_XY pp1[2] = { p1, gp_XY( p1.X() - grad.Y(),
+ p1.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp1, 0, POS_LEFT|POS_RIGHT )._name == POS_LEFT )
+ return PointPos( POS_VERTEX, index0 );
+
+ // check if the point2D is at "vertex 2" zone
+ gp_XY pp2[2] = { p2, gp_XY( p2.X() - grad.Y(),
+ p2.Y() + grad.X() ) };
+ if ( getPointPosition( point2D, pp2, 0, POS_LEFT|POS_RIGHT )._name == POS_RIGHT )
+ return PointPos( POS_VERTEX, index0 + 1);
+ }
+ double edgeEquation =
+ ( point2D.X() - p1.X() ) * grad.Y() - ( point2D.Y() - p1.Y() ) * grad.X();
+ return PointPos( edgeEquation < 0 ? POS_LEFT : POS_RIGHT, index0 );
+ }
+}
+
+//=======================================================================
+/*!
+ * \brief Return minimal distance from a point to a face
+ *
+ * Currently we ignore non-planarity and 2nd order of face
+ */
+//=======================================================================
+
+double SMESH_MeshEditor::GetDistance( const SMDS_MeshFace* face,
+ const gp_Pnt& point )
+{
+ double badDistance = -1;
+ if ( !face ) return badDistance;
+
+ // coordinates of nodes (medium nodes, if any, ignored)
+ typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator;
+ vector<gp_XYZ> xyz( TXyzIterator( face->nodesIterator()), TXyzIterator() );
+ xyz.resize( face->NbCornerNodes()+1 );
+
+ // transformation to get xyz[0] lies on the origin, xyz[1] lies on the Z axis,
+ // and xyz[2] lies in the XZ plane. This is to pass to 2D space on XZ plane.
+ gp_Trsf trsf;
+ gp_Vec OZ ( xyz[0], xyz[1] );
+ gp_Vec OX ( xyz[0], xyz[2] );
+ if ( OZ.Magnitude() < std::numeric_limits<double>::min() )
+ {
+ if ( xyz.size() < 4 ) return badDistance;
+ OZ = gp_Vec ( xyz[0], xyz[2] );
+ OX = gp_Vec ( xyz[0], xyz[3] );
+ }
+ gp_Ax3 tgtCS;
+ try {
+ tgtCS = gp_Ax3( xyz[0], OZ, OX );
+ }
+ catch ( Standard_Failure ) {
+ return badDistance;
+ }
+ trsf.SetTransformation( tgtCS );
+
+ // move all the nodes to 2D
+ vector<gp_XY> xy( xyz.size() );
+ for ( size_t i = 0;i < xyz.size()-1; ++i )
+ {
+ gp_XYZ p3d = xyz[i];
+ trsf.Transforms( p3d );
+ xy[i].SetCoord( p3d.X(), p3d.Z() );
+ }
+ xyz.back() = xyz.front();
+ xy.back() = xy.front();
+
+ // // move the point in 2D
+ gp_XYZ tmpPnt = point.XYZ();
+ trsf.Transforms( tmpPnt );
+ gp_XY point2D( tmpPnt.X(), tmpPnt.Z() );
+
+ // loop on segments of the face to analyze point position ralative to the face
+ set< PointPos > pntPosSet;
+ for ( size_t i = 1; i < xy.size(); ++i )
+ {
+ PointPos pos = getPointPosition( point2D, &xy[0], i-1 );
+ pntPosSet.insert( pos );
+ }
+
+ // compute distance
+ PointPos pos = *pntPosSet.begin();
+ // cout << "Face " << face->GetID() << " DIST: ";
+ switch ( pos._name )
+ {
+ case POS_LEFT: {
+ // point is most close to a segment
+ gp_Vec p0p1( point, xyz[ pos._index ] );
+ gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector
+ p1p2.Normalize();
+ double projDist = p0p1 * p1p2; // distance projected to the segment
+ gp_Vec projVec = p1p2 * projDist;
+ gp_Vec distVec = p0p1 - projVec;
+ // cout << distVec.Magnitude() << ", SEG " << face->GetNode(pos._index)->GetID()
+ // << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ case POS_RIGHT: {
+ // point is inside the face
+ double distToFacePlane = tmpPnt.Y();
+ // cout << distToFacePlane << ", INSIDE " << endl;
+ return Abs( distToFacePlane );
+ }
+ case POS_VERTEX: {
+ // point is most close to a node
+ gp_Vec distVec( point, xyz[ pos._index ]);
+ // cout << distVec.Magnitude() << " VERTEX " << face->GetNode(pos._index)->GetID() << endl;
+ return distVec.Magnitude();
+ }
+ }
+ return badDistance;
+}
+
//=======================================================================
//function : SimplifyFace
//purpose :
SMESHDS_Mesh* aMesh = GetMeshDS();
// TODO algoritm not OK with vtkUnstructuredGrid: 2 meshes can't share nodes
//SMDS_Mesh aTmpFacesMesh; // try to use the same mesh
- set<const SMDS_MeshElement*> faceSet1, faceSet2;
+ TIDSortedElemSet faceSet1, faceSet2;
set<const SMDS_MeshElement*> volSet1, volSet2;
set<const SMDS_MeshNode*> nodeSet1, nodeSet2;
- set<const SMDS_MeshElement*> * faceSetPtr[] = { &faceSet1, &faceSet2 };
- set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
+ TIDSortedElemSet * faceSetPtr[] = { &faceSet1, &faceSet2 };
+ set<const SMDS_MeshElement*> * volSetPtr[] = { &volSet1, &volSet2 };
set<const SMDS_MeshNode*> * nodeSetPtr[] = { &nodeSet1, &nodeSet2 };
- TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
+ TIDSortedElemSet * elemSetPtr[] = { &theSide1, &theSide2 };
int iSide, iFace, iNode;
list<const SMDS_MeshElement* > tempFaceList;
for ( iSide = 0; iSide < 2; iSide++ ) {
set<const SMDS_MeshNode*> * nodeSet = nodeSetPtr[ iSide ];
- TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
- set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
+ TIDSortedElemSet * elemSet = elemSetPtr[ iSide ];
+ TIDSortedElemSet * faceSet = faceSetPtr[ iSide ];
set<const SMDS_MeshElement*> * volSet = volSetPtr [ iSide ];
set<const SMDS_MeshElement*>::iterator vIt;
TIDSortedElemSet::iterator eIt;
bool isNewFace = setOfFaceNodeSet.insert( faceNodeSet ).second;
if ( isNewFace ) {
// no such a face is given but it still can exist, check it
- if ( nbNodes == 3 ) {
- aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2] );
- }
- else if ( nbNodes == 4 ) {
- aFreeFace = aMesh->FindFace( fNodes[0],fNodes[1],fNodes[2],fNodes[3] );
- }
- else {
- vector<const SMDS_MeshNode *> poly_nodes ( fNodes, & fNodes[nbNodes]);
- aFreeFace = aMesh->FindFace(poly_nodes);
- }
+ vector<const SMDS_MeshNode *> nodes ( fNodes, fNodes + nbNodes);
+ aFreeFace = aMesh->FindElement( nodes, SMDSAbs_Face, /*noMedium=*/false );
}
if ( !aFreeFace ) {
// create a temporary face
//aFreeFace = aTmpFacesMesh.AddPolygonalFace(poly_nodes);
aFreeFace = aMesh->AddPolygonalFace(poly_nodes);
}
+ if ( aFreeFace )
+ tempFaceList.push_back( aFreeFace );
}
- if ( aFreeFace ) {
+
+ if ( aFreeFace )
freeFaceList.push_back( aFreeFace );
- tempFaceList.push_back( aFreeFace );
- }
} // loop on faces of a volume
- // choose one of several free faces
- // --------------------------------------
+ // choose one of several free faces of a volume
+ // --------------------------------------------
if ( freeFaceList.size() > 1 ) {
// choose a face having max nb of nodes shared by other elems of a side
- int maxNbNodes = -1/*, nbExcludedFaces = 0*/;
+ int maxNbNodes = -1;
list<const SMDS_MeshElement* >::iterator fIt = freeFaceList.begin();
while ( fIt != freeFaceList.end() ) { // loop on free faces
int nbSharedNodes = 0;
SMDS_ElemIteratorPtr invElemIt = n->GetInverseElementIterator();
while ( invElemIt->more() ) {
const SMDS_MeshElement* e = invElemIt->next();
- if ( faceSet->find( e ) != faceSet->end() )
- nbSharedNodes++;
- if ( elemSet->find( e ) != elemSet->end() )
- nbSharedNodes++;
+ nbSharedNodes += faceSet->count( e );
+ nbSharedNodes += elemSet->count( e );
}
}
- if ( nbSharedNodes >= maxNbNodes ) {
+ if ( nbSharedNodes > maxNbNodes ) {
maxNbNodes = nbSharedNodes;
+ freeFaceList.erase( freeFaceList.begin(), fIt++ );
+ }
+ else if ( nbSharedNodes == maxNbNodes ) {
fIt++;
}
- else
+ else {
freeFaceList.erase( fIt++ ); // here fIt++ occurs before erase
+ }
}
if ( freeFaceList.size() > 1 )
{
TNodeNodeMap nReplaceMap; // bind a node to remove to a node to put instead
if ( theFirstNode1 != theFirstNode2 )
- nReplaceMap.insert( TNodeNodeMap::value_type( theFirstNode1, theFirstNode2 ));
+ nReplaceMap.insert( make_pair( theFirstNode1, theFirstNode2 ));
if ( theSecondNode1 != theSecondNode2 )
- nReplaceMap.insert( TNodeNodeMap::value_type( theSecondNode1, theSecondNode2 ));
+ nReplaceMap.insert( make_pair( theSecondNode1, theSecondNode2 ));
LinkID_Gen aLinkID_Gen( GetMeshDS() );
set< long > linkIdSet; // links to process
// loop on links in linkList; find faces by links and append links
// of the found faces to linkList
list< NLink >::iterator linkIt[] = { linkList[0].begin(), linkList[1].begin() } ;
- for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ ) {
+ for ( ; linkIt[0] != linkList[0].end(); linkIt[0]++, linkIt[1]++ )
+ {
NLink link[] = { *linkIt[0], *linkIt[1] };
long linkID = aLinkID_Gen.GetLinkID( link[0].first, link[0].second );
- if ( linkIdSet.find( linkID ) == linkIdSet.end() )
+ if ( !linkIdSet.count( linkID ) )
continue;
// by links, find faces in the face sets,
// ---------------------------------------------------------------
const SMDS_MeshElement* face[] = { 0, 0 };
- //const SMDS_MeshNode* faceNodes[ 2 ][ 5 ];
- vector<const SMDS_MeshNode*> fnodes1(9);
- vector<const SMDS_MeshNode*> fnodes2(9);
- //const SMDS_MeshNode* notLinkNodes[ 2 ][ 2 ] = {{ 0, 0 },{ 0, 0 }} ;
- vector<const SMDS_MeshNode*> notLinkNodes1(6);
- vector<const SMDS_MeshNode*> notLinkNodes2(6);
+ vector<const SMDS_MeshNode*> fnodes[2];
int iLinkNode[2][2];
+ TIDSortedElemSet avoidSet;
for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
const SMDS_MeshNode* n1 = link[iSide].first;
const SMDS_MeshNode* n2 = link[iSide].second;
- set<const SMDS_MeshElement*> * faceSet = faceSetPtr[ iSide ];
- set< const SMDS_MeshElement* > fMap;
- for ( int i = 0; i < 2; i++ ) { // loop on 2 nodes of a link
- const SMDS_MeshNode* n = i ? n1 : n2; // a node of a link
- SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face);
- while ( fIt->more() ) { // loop on faces sharing a node
- const SMDS_MeshElement* f = fIt->next();
- if (faceSet->find( f ) != faceSet->end() && // f is in face set
- ! fMap.insert( f ).second ) // f encounters twice
- {
- if ( face[ iSide ] ) {
- MESSAGE( "2 faces per link " );
- aResult = iSide ? SEW_BAD_SIDE2_NODES : SEW_BAD_SIDE1_NODES;
- break;
- }
- face[ iSide ] = f;
- faceSet->erase( f );
- // get face nodes and find ones of a link
- iNode = 0;
- int nbl = -1;
- if(f->IsPoly()) {
- if(iSide==0) {
- fnodes1.resize(f->NbNodes()+1);
- notLinkNodes1.resize(f->NbNodes()-2);
- }
- else {
- fnodes2.resize(f->NbNodes()+1);
- notLinkNodes2.resize(f->NbNodes()-2);
- }
- }
- if(!f->IsQuadratic()) {
- SMDS_ElemIteratorPtr nIt = f->nodesIterator();
- while ( nIt->more() ) {
- const SMDS_MeshNode* n =
- static_cast<const SMDS_MeshNode*>( nIt->next() );
- if ( n == n1 ) {
- iLinkNode[ iSide ][ 0 ] = iNode;
- }
- else if ( n == n2 ) {
- iLinkNode[ iSide ][ 1 ] = iNode;
- }
- //else if ( notLinkNodes[ iSide ][ 0 ] )
- // notLinkNodes[ iSide ][ 1 ] = n;
- //else
- // notLinkNodes[ iSide ][ 0 ] = n;
- else {
- nbl++;
- if(iSide==0)
- notLinkNodes1[nbl] = n;
- //notLinkNodes1.push_back(n);
- else
- notLinkNodes2[nbl] = n;
- //notLinkNodes2.push_back(n);
- }
- //faceNodes[ iSide ][ iNode++ ] = n;
- if(iSide==0) {
- fnodes1[iNode++] = n;
- }
- else {
- fnodes2[iNode++] = n;
- }
- }
- }
- else { // f->IsQuadratic()
- const SMDS_VtkFace* F =
- dynamic_cast<const SMDS_VtkFace*>(f);
- if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
- // use special nodes iterator
- SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
- while ( anIter->more() ) {
- const SMDS_MeshNode* n =
- static_cast<const SMDS_MeshNode*>( anIter->next() );
- if ( n == n1 ) {
- iLinkNode[ iSide ][ 0 ] = iNode;
- }
- else if ( n == n2 ) {
- iLinkNode[ iSide ][ 1 ] = iNode;
- }
- else {
- nbl++;
- if(iSide==0) {
- notLinkNodes1[nbl] = n;
- }
- else {
- notLinkNodes2[nbl] = n;
- }
- }
- if(iSide==0) {
- fnodes1[iNode++] = n;
- }
- else {
- fnodes2[iNode++] = n;
- }
- }
- }
- //faceNodes[ iSide ][ iNode ] = faceNodes[ iSide ][ 0 ];
- if(iSide==0) {
- fnodes1[iNode] = fnodes1[0];
- }
- else {
- fnodes2[iNode] = fnodes1[0];
- }
- }
+ //cout << "Side " << iSide << " ";
+ //cout << "L( " << n1->GetID() << ", " << n2->GetID() << " ) " << endl;
+ // find a face by two link nodes
+ face[ iSide ] = FindFaceInSet( n1, n2, *faceSetPtr[ iSide ], avoidSet,
+ &iLinkNode[iSide][0], &iLinkNode[iSide][1] );
+ if ( face[ iSide ])
+ {
+ //cout << " F " << face[ iSide]->GetID() <<endl;
+ faceSetPtr[ iSide ]->erase( face[ iSide ]);
+ // put face nodes to fnodes
+ if ( face[ iSide ]->IsQuadratic() )
+ {
+ // use interlaced nodes iterator
+ const SMDS_VtkFace* F = dynamic_cast<const SMDS_VtkFace*>( face[ iSide ]);
+ if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace"));
+ SMDS_ElemIteratorPtr nIter = F->interlacedNodesElemIterator();
+ while ( nIter->more() )
+ fnodes[ iSide ].push_back( cast2Node( nIter->next() ));
}
+ else
+ {
+ fnodes[ iSide ].assign( face[ iSide ]->begin_nodes(),
+ face[ iSide ]->end_nodes() );
+ }
+ fnodes[ iSide ].push_back( fnodes[ iSide ].front());
}
}
else {
aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
}
- break; // do not return because it s necessary to remove tmp faces
+ break; // do not return because it's necessary to remove tmp faces
}
// set nodes to merge
// -------------------
if ( face[0] && face[1] ) {
- int nbNodes = face[0]->NbNodes();
+ const int nbNodes = face[0]->NbNodes();
if ( nbNodes != face[1]->NbNodes() ) {
MESSAGE("Diff nb of face nodes");
aResult = SEW_TOPO_DIFF_SETS_OF_ELEMENTS;
break; // do not return because it s necessary to remove tmp faces
}
- bool reverse[] = { false, false }; // order of notLinkNodes of quadrangle
- if ( nbNodes == 3 ) {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[0], notLinkNodes2[0] ));
+ bool reverse[] = { false, false }; // order of nodes in the link
+ for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
+ // analyse link orientation in faces
+ int i1 = iLinkNode[ iSide ][ 0 ];
+ int i2 = iLinkNode[ iSide ][ 1 ];
+ reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
}
- else {
- for ( iSide = 0; iSide < 2; iSide++ ) { // loop on 2 sides
- // analyse link orientation in faces
- int i1 = iLinkNode[ iSide ][ 0 ];
- int i2 = iLinkNode[ iSide ][ 1 ];
- reverse[ iSide ] = Abs( i1 - i2 ) == 1 ? i1 > i2 : i2 > i1;
- // if notLinkNodes are the first and the last ones, then
- // their order does not correspond to the link orientation
- if (( i1 == 1 && i2 == 2 ) ||
- ( i1 == 2 && i2 == 1 ))
- reverse[ iSide ] = !reverse[ iSide ];
- }
- if ( reverse[0] == reverse[1] ) {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][0] ));
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][1], notLinkNodes[1][1] ));
- for(int nn=0; nn<nbNodes-2; nn++) {
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[nn], notLinkNodes2[nn] ));
- }
- }
- else {
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][0], notLinkNodes[1][1] ));
- //nReplaceMap.insert( TNodeNodeMap::value_type
- // ( notLinkNodes[0][1], notLinkNodes[1][0] ));
- for(int nn=0; nn<nbNodes-2; nn++) {
- nReplaceMap.insert( TNodeNodeMap::value_type
- ( notLinkNodes1[nn], notLinkNodes2[nbNodes-3-nn] ));
- }
- }
+ int di1 = reverse[0] ? -1 : +1, i1 = iLinkNode[0][1] + di1;
+ int di2 = reverse[1] ? -1 : +1, i2 = iLinkNode[1][1] + di2;
+ for ( int i = nbNodes - 2; i > 0; --i, i1 += di1, i2 += di2 )
+ {
+ nReplaceMap.insert ( make_pair ( fnodes[0][ ( i1 + nbNodes ) % nbNodes ],
+ fnodes[1][ ( i2 + nbNodes ) % nbNodes ]));
}
// add other links of the faces to linkList
// -----------------------------------------
- //const SMDS_MeshNode** nodes = faceNodes[ 0 ];
for ( iNode = 0; iNode < nbNodes; iNode++ ) {
- //linkID = aLinkID_Gen.GetLinkID( nodes[iNode], nodes[iNode+1] );
- linkID = aLinkID_Gen.GetLinkID( fnodes1[iNode], fnodes1[iNode+1] );
+ linkID = aLinkID_Gen.GetLinkID( fnodes[0][iNode], fnodes[0][iNode+1] );
pair< set<long>::iterator, bool > iter_isnew = linkIdSet.insert( linkID );
if ( !iter_isnew.second ) { // already in a set: no need to process
linkIdSet.erase( iter_isnew.first );
}
else // new in set == encountered for the first time: add
{
- //const SMDS_MeshNode* n1 = nodes[ iNode ];
- //const SMDS_MeshNode* n2 = nodes[ iNode + 1];
- const SMDS_MeshNode* n1 = fnodes1[ iNode ];
- const SMDS_MeshNode* n2 = fnodes1[ iNode + 1];
+ const SMDS_MeshNode* n1 = fnodes[0][ iNode ];
+ const SMDS_MeshNode* n2 = fnodes[0][ iNode + 1];
linkList[0].push_back ( NLink( n1, n2 ));
linkList[1].push_back ( NLink( nReplaceMap[n1], nReplaceMap[n2] ));
}
}
} // 2 faces found
+
+ if ( faceSetPtr[0]->empty() || faceSetPtr[1]->empty() )
+ break;
+
} // loop on link lists
if ( aResult == SEW_OK &&
- ( linkIt[0] != linkList[0].end() ||
+ ( //linkIt[0] != linkList[0].end() ||
!faceSetPtr[0]->empty() || !faceSetPtr[1]->empty() )) {
MESSAGE( (linkIt[0] != linkList[0].end()) <<" "<< (faceSetPtr[0]->empty()) <<
" " << (faceSetPtr[1]->empty()));
// 3. Replace nodes in elements of the side 1 and remove replaced nodes
// ====================================================================
- // delete temporary faces: they are in reverseElements of actual nodes
+ // delete temporary faces
// SMDS_FaceIteratorPtr tmpFaceIt = aTmpFacesMesh.facesIterator();
// while ( tmpFaceIt->more() )
// aTmpFacesMesh.RemoveElement( tmpFaceIt->next() );
-// list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
-// for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
-// aMesh->RemoveElement(*tmpFaceIt);
+ list<const SMDS_MeshElement* >::iterator tmpFaceIt = tempFaceList.begin();
+ for (; tmpFaceIt !=tempFaceList.end(); ++tmpFaceIt)
+ aMesh->RemoveElement(*tmpFaceIt);
if ( aResult != SEW_OK)
return aResult;
