X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=1d5698d249bc67e59ed207f38ab86ebc735811ca;hp=b9a15a21ee12671f231aac27d9f3d98da6332e52;hb=3f3beb8b958607b3ac214b8234b5ae5d99facc13;hpb=9c2ce3dffb1f26145f6ab49ebe9c9bb61254725b diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index b9a15a21e..1d5698d24 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -42,6 +42,7 @@ #include "SMESH_Algo.hxx" #include "SMESH_ControlsDef.hxx" #include "SMESH_Group.hxx" +#include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_OctreeNode.hxx" #include "SMESH_subMesh.hxx" @@ -58,15 +59,10 @@ #include #include #include -#include #include -#include #include #include -#include #include -#include -#include #include #include #include @@ -107,10 +103,16 @@ using namespace std; using namespace SMESH::Controls; -typedef map > TElemOfNodeListMap; -typedef map > TElemOfElemListMap; - -typedef SMDS_SetIterator< SMDS_pElement, TIDSortedElemSet::const_iterator> TSetIterator; +namespace +{ + template < class ELEM_SET > + SMDS_ElemIteratorPtr elemSetIterator( const ELEM_SET& elements ) + { + typedef SMDS_SetIterator + < SMDS_pElement, typename ELEM_SET::const_iterator> TSetIterator; + return SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() )); + } +} //======================================================================= //function : SMESH_MeshEditor @@ -169,6 +171,12 @@ SMESH_MeshEditor::AddElement(const vector & node, else e = mesh->AddFace (node[0], node[1], node[2], node[3], node[4], node[5] ); } + else if (nbnode == 7) { + if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], ID); + else e = mesh->AddFace (node[0], node[1], node[2], node[3], + node[4], node[5], node[6] ); + } else if (nbnode == 8) { if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], node[4], node[5], node[6], node[7], ID); @@ -417,12 +425,21 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs, void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements, TIDSortedElemSet& all0DElems ) { - typedef SMDS_SetIterator TSetIterator; SMDS_ElemIteratorPtr elemIt; + vector< const SMDS_MeshElement* > allNodes; if ( elements.empty() ) + { + allNodes.reserve( GetMeshDS()->NbNodes() ); elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node ); + while ( elemIt->more() ) + allNodes.push_back( elemIt->next() ); + + elemIt = elemSetIterator( allNodes ); + } else - elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() )); + { + elemIt = elemSetIterator( elements ); + } while ( elemIt->more() ) { @@ -530,12 +547,12 @@ bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode* node, } //======================================================================= -//function : ShiftNodesQuadTria -//purpose : auxilary -// Shift nodes in the array corresponded to quadratic triangle +//function : shiftNodesQuadTria +//purpose : Shift nodes in the array corresponded to quadratic triangle // example: (0,1,2,3,4,5) -> (1,2,0,4,5,3) //======================================================================= -static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[]) + +static void shiftNodesQuadTria(vector< const SMDS_MeshNode* >& aNodes) { const SMDS_MeshNode* nd1 = aNodes[0]; aNodes[0] = aNodes[1]; @@ -548,53 +565,42 @@ static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[]) } //======================================================================= -//function : edgeConnectivity -//purpose : auxilary -// return number of the edges connected with the theNode. +//function : nbEdgeConnectivity +//purpose : return number of the edges connected with the theNode. // if theEdges has connections with the other type of the // elements, return -1 //======================================================================= + static int nbEdgeConnectivity(const SMDS_MeshNode* theNode) { - SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(); - int nb=0; - while(elemIt->more()) { - elemIt->next(); - nb++; - } - return nb; + // SMDS_ElemIteratorPtr elemIt = theNode->GetInverseElementIterator(); + // int nb=0; + // while(elemIt->more()) { + // elemIt->next(); + // nb++; + // } + // return nb; + return theNode->NbInverseElements(); } - //======================================================================= -//function : GetNodesFromTwoTria -//purpose : auxilary -// Shift nodes in the array corresponded to quadratic triangle -// example: (0,1,2,3,4,5) -> (1,2,0,4,5,3) +//function : getNodesFromTwoTria +//purpose : //======================================================================= -static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1, + +static bool getNodesFromTwoTria(const SMDS_MeshElement * theTria1, const SMDS_MeshElement * theTria2, - const SMDS_MeshNode* N1[], - const SMDS_MeshNode* N2[]) + vector< const SMDS_MeshNode*>& N1, + vector< const SMDS_MeshNode*>& N2) { - SMDS_ElemIteratorPtr it = theTria1->nodesIterator(); - int i=0; - while(i<6) { - N1[i] = static_cast( it->next() ); - i++; - } - if(it->more()) return false; - it = theTria2->nodesIterator(); - i=0; - while(i<6) { - N2[i] = static_cast( it->next() ); - i++; - } - if(it->more()) return false; + N1.assign( theTria1->begin_nodes(), theTria1->end_nodes() ); + if ( N1.size() < 6 ) return false; + N2.assign( theTria2->begin_nodes(), theTria2->end_nodes() ); + if ( N2.size() < 6 ) return false; int sames[3] = {-1,-1,-1}; int nbsames = 0; - int j; + int i, j; for(i=0; i<3; i++) { for(j=0; j<3; j++) { if(N1[i]==N2[j]) { @@ -606,18 +612,18 @@ static bool GetNodesFromTwoTria(const SMDS_MeshElement * theTria1, } if(nbsames!=2) return false; if(sames[0]>-1) { - ShiftNodesQuadTria(N1); + shiftNodesQuadTria(N1); if(sames[1]>-1) { - ShiftNodesQuadTria(N1); + shiftNodesQuadTria(N1); } } i = sames[0] + sames[1] + sames[2]; for(; i<2; i++) { - ShiftNodesQuadTria(N2); + shiftNodesQuadTria(N2); } - // now we receive following N1 and N2 (using numeration as above image) + // now we receive following N1 and N2 (using numeration as in the image below) // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6) - // i.e. first nodes from both arrays determ new diagonal + // i.e. first nodes from both arrays form a new diagonal return true; } @@ -706,9 +712,11 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1, } // end if(F1 && F2) // check case of quadratic faces - if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle) + if (theTria1->GetEntityType() != SMDSEntity_Quad_Triangle && + theTria1->GetEntityType() != SMDSEntity_BiQuad_Triangle) return false; - if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle) + if (theTria2->GetEntityType() != SMDSEntity_Quad_Triangle&& + theTria2->GetEntityType() != SMDSEntity_BiQuad_Triangle) return false; // 5 @@ -721,32 +729,61 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1, // 4 +--+--+ 3 // 8 - const SMDS_MeshNode* N1 [6]; - const SMDS_MeshNode* N2 [6]; - if(!GetNodesFromTwoTria(theTria1,theTria2,N1,N2)) + vector< const SMDS_MeshNode* > N1; + vector< const SMDS_MeshNode* > N2; + if(!getNodesFromTwoTria(theTria1,theTria2,N1,N2)) return false; // now we receive following N1 and N2 (using numeration as above image) // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6) // i.e. first nodes from both arrays determ new diagonal - const SMDS_MeshNode* N1new [6]; - const SMDS_MeshNode* N2new [6]; - N1new[0] = N1[0]; - N1new[1] = N2[0]; - N1new[2] = N2[1]; - N1new[3] = N1[4]; - N1new[4] = N2[3]; - N1new[5] = N1[5]; - N2new[0] = N1[0]; - N2new[1] = N1[1]; - N2new[2] = N2[0]; - N2new[3] = N1[3]; - N2new[4] = N2[5]; - N2new[5] = N1[4]; - // replaces nodes in faces - GetMeshDS()->ChangeElementNodes( theTria1, N1new, 6 ); - GetMeshDS()->ChangeElementNodes( theTria2, N2new, 6 ); - + vector< const SMDS_MeshNode*> N1new( N1.size() ); + vector< const SMDS_MeshNode*> N2new( N2.size() ); + N1new.back() = N1.back(); // central node of biquadratic + N2new.back() = N2.back(); + N1new[0] = N1[0]; N2new[0] = N1[0]; + N1new[1] = N2[0]; N2new[1] = N1[1]; + N1new[2] = N2[1]; N2new[2] = N2[0]; + N1new[3] = N1[4]; N2new[3] = N1[3]; + N1new[4] = N2[3]; N2new[4] = N2[5]; + N1new[5] = N1[5]; N2new[5] = N1[4]; + // change nodes in faces + GetMeshDS()->ChangeElementNodes( theTria1, &N1new[0], N1new.size() ); + GetMeshDS()->ChangeElementNodes( theTria2, &N2new[0], N2new.size() ); + + // move the central node of biquadratic triangle + SMESH_MesherHelper helper( *GetMesh() ); + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + const SMDS_MeshElement* tria = is2nd ? theTria2 : theTria1; + vector< const SMDS_MeshNode*>& nodes = is2nd ? N2new : N1new; + if ( nodes.size() < 7 ) + continue; + helper.SetSubShape( tria->getshapeId() ); + const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() ); + gp_Pnt xyz; + if ( F.IsNull() ) + { + xyz = ( SMESH_TNodeXYZ( nodes[3] ) + + SMESH_TNodeXYZ( nodes[4] ) + + SMESH_TNodeXYZ( nodes[5] )) / 3.; + } + else + { + bool checkUV; + gp_XY uv = ( helper.GetNodeUV( F, nodes[3], nodes[2], &checkUV ) + + helper.GetNodeUV( F, nodes[4], nodes[0], &checkUV ) + + helper.GetNodeUV( F, nodes[5], nodes[1], &checkUV )) / 3.; + TopLoc_Location loc; + Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc); + xyz = S->Value( uv.X(), uv.Y() ); + xyz.Transform( loc ); + if ( nodes[6]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE && // set UV + nodes[6]->getshapeId() > 0 ) + GetMeshDS()->SetNodeOnFace( nodes[6], nodes[6]->getshapeId(), uv.X(), uv.Y() ); + } + GetMeshDS()->MoveNode( nodes[6], xyz.X(), xyz.Y(), xyz.Z() ); + } return true; } @@ -768,30 +805,28 @@ static bool findTriangles(const SMDS_MeshNode * theNode1, SMDS_ElemIteratorPtr it = theNode1->GetInverseElementIterator(SMDSAbs_Face); while (it->more()) { const SMDS_MeshElement* elem = it->next(); - if ( elem->NbNodes() == 3 ) + if ( elem->NbCornerNodes() == 3 ) emap.insert( elem ); } it = theNode2->GetInverseElementIterator(SMDSAbs_Face); while (it->more()) { const SMDS_MeshElement* elem = it->next(); - if ( emap.find( elem ) != emap.end() ) { - if ( theTria1 ) { - // theTria1 must be element with minimum ID - if( theTria1->GetID() < elem->GetID() ) { - theTria2 = elem; - } - else { - theTria2 = theTria1; - theTria1 = elem; - } - break; - } - else { + if ( emap.count( elem )) { + if ( !theTria1 ) + { theTria1 = elem; } + else + { + theTria2 = elem; + // theTria1 must be element with minimum ID + if ( theTria2->GetID() < theTria1->GetID() ) + std::swap( theTria2, theTria1 ); + return true; + } } } - return ( theTria1 && theTria2 ); + return false; } //======================================================================= @@ -980,9 +1015,9 @@ bool SMESH_MeshEditor::DeleteDiag (const SMDS_MeshNode * theNode1, // 4 +--+--+ 3 // 8 - const SMDS_MeshNode* N1 [6]; - const SMDS_MeshNode* N2 [6]; - if(!GetNodesFromTwoTria(tr1,tr2,N1,N2)) + vector< const SMDS_MeshNode* > N1; + vector< const SMDS_MeshNode* > N2; + if(!getNodesFromTwoTria(tr1,tr2,N1,N2)) return false; // now we receive following N1 and N2 (using numeration as above image) // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6) @@ -1034,82 +1069,49 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem) if ( !it || !it->more() ) return false; - switch ( theElem->GetType() ) { + const SMDSAbs_ElementType type = theElem->GetType(); + if ( type < SMDSAbs_Edge || type > SMDSAbs_Volume ) + return false; - case SMDSAbs_Edge: - case SMDSAbs_Face: { - if(!theElem->IsQuadratic()) { - int i = theElem->NbNodes(); - vector aNodes( i ); - while ( it->more() ) - aNodes[ --i ]= static_cast( it->next() ); - return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], theElem->NbNodes() ); - } - else { - // quadratic elements - if(theElem->GetType()==SMDSAbs_Edge) { - vector aNodes(3); - aNodes[1]= static_cast( it->next() ); - aNodes[0]= static_cast( it->next() ); - aNodes[2]= static_cast( it->next() ); - return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], 3 ); - } - else { - int nbn = theElem->NbNodes(); - vector aNodes(nbn); - aNodes[0]= static_cast( it->next() ); - int i=1; - for(; i( it->next() ); - } - for(i=0; i( it->next() ); - } - return GetMeshDS()->ChangeElementNodes( theElem, &aNodes[0], nbn ); - } + const SMDSAbs_EntityType geomType = theElem->GetEntityType(); + if ( geomType == SMDSEntity_Polyhedra ) // polyhedron + { + const SMDS_VtkVolume* aPolyedre = + dynamic_cast( theElem ); + if (!aPolyedre) { + MESSAGE("Warning: bad volumic element"); + return false; } - } - case SMDSAbs_Volume: { - if (theElem->IsPoly()) { - // TODO reorient vtk polyhedron - MESSAGE("reorient vtk polyhedron ?"); - const SMDS_VtkVolume* aPolyedre = - dynamic_cast( theElem ); - if (!aPolyedre) { - MESSAGE("Warning: bad volumic element"); - return false; - } - - int nbFaces = aPolyedre->NbFaces(); - vector poly_nodes; - vector quantities (nbFaces); + const int nbFaces = aPolyedre->NbFaces(); + vector poly_nodes; + vector quantities (nbFaces); - // reverse each face of the polyedre - for (int iface = 1; iface <= nbFaces; iface++) { - int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface); - quantities[iface - 1] = nbFaceNodes; + // reverse each face of the polyedre + for (int iface = 1; iface <= nbFaces; iface++) { + int inode, nbFaceNodes = aPolyedre->NbFaceNodes(iface); + quantities[iface - 1] = nbFaceNodes; - for (inode = nbFaceNodes; inode >= 1; inode--) { - const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode); - poly_nodes.push_back(curNode); - } + for (inode = nbFaceNodes; inode >= 1; inode--) { + const SMDS_MeshNode* curNode = aPolyedre->GetFaceNode(iface, inode); + poly_nodes.push_back(curNode); } - - return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities ); - - } - else { - SMDS_VolumeTool vTool; - if ( !vTool.Set( theElem )) - return false; - vTool.Inverse(); - MESSAGE("ChangeElementNodes reorient: check vTool.Inverse"); - return GetMeshDS()->ChangeElementNodes( theElem, vTool.GetNodes(), vTool.NbNodes() ); } + return GetMeshDS()->ChangePolyhedronNodes( theElem, poly_nodes, quantities ); } - default:; + else // other elements + { + vector nodes( theElem->begin_nodes(), theElem->end_nodes() ); + const std::vector& interlace = SMDS_MeshCell::reverseSmdsOrder( geomType ); + if ( interlace.empty() ) + { + std::reverse( nodes.begin(), nodes.end() ); // polygon + } + else if ( interlace.size() > 1 ) + { + SMDS_MeshCell::applyInterlace( interlace, nodes ); + } + return GetMeshDS()->ChangeElementNodes( theElem, &nodes[0], nodes.size() ); } - return false; } @@ -1140,7 +1142,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, // orient theFace according to theDirection gp_XYZ normal; - SMESH_Algo::FaceNormal( theFace, normal, /*normalized=*/false ); + SMESH_MeshAlgos::FaceNormal( theFace, normal, /*normalized=*/false ); if ( normal * theDirection.XYZ() < 0 ) nbReori += Reorient( theFace ); @@ -1189,8 +1191,9 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, { facesNearLink.clear(); nodeIndsOfFace.clear(); - while (( otherFace = FindFaceInSet( link.first, link.second, - theFaces, avoidSet, &nodeInd1, &nodeInd2 ))) + while (( otherFace = SMESH_MeshAlgos::FindFaceInSet( link.first, link.second, + theFaces, avoidSet, + &nodeInd1, &nodeInd2 ))) if ( otherFace != theFace) { facesNearLink.push_back( otherFace ); @@ -1203,10 +1206,10 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, // select a face most co-directed with theFace, // other faces won't be visited this time gp_XYZ NF, NOF; - SMESH_Algo::FaceNormal( theFace, NF, /*normalized=*/false ); + SMESH_MeshAlgos::FaceNormal( theFace, NF, /*normalized=*/false ); double proj, maxProj = -1; for ( size_t i = 0; i < facesNearLink.size(); ++i ) { - SMESH_Algo::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false ); + SMESH_MeshAlgos::FaceNormal( facesNearLink[i], NOF, /*normalized=*/false ); if (( proj = Abs( NF * NOF )) > maxProj ) { maxProj = proj; otherFace = facesNearLink[i]; @@ -1269,8 +1272,6 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, myLastCreatedElems.Clear(); myLastCreatedNodes.Clear(); - MESSAGE( "::QuadToTri()" ); - if ( !theCrit.get() ) return false; @@ -1280,7 +1281,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, SMESH_MesherHelper helper( *GetMesh() ); TIDSortedElemSet::iterator itElem; - for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) { + for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) + { const SMDS_MeshElement* elem = *itElem; if ( !elem || elem->GetType() != SMDSAbs_Face ) continue; @@ -1300,13 +1302,12 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] ); aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit ); - int aShapeId = FindShape( elem ); + const int aShapeId = FindShape( elem ); const SMDS_MeshElement* newElem1 = 0; const SMDS_MeshElement* newElem2 = 0; - if( !elem->IsQuadratic() ) { - - // split liner quadrangle + if ( !elem->IsQuadratic() ) // split liner quadrangle + { // for MaxElementLength2D functor we return minimum diagonal for splitting, // because aBadRate1=2*len(diagonal 1-3); aBadRate2=2*len(diagonal 2-4) if ( aBadRate1 <= aBadRate2 ) { @@ -1320,70 +1321,28 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, newElem2 = aMesh->AddFace( aNodes[3], aNodes[1], aNodes[2] ); } } - else { - - // split quadratic quadrangle + else // split quadratic quadrangle + { + helper.SetIsQuadratic( true ); + helper.SetIsBiQuadratic( aNodes.size() == 9 ); - // get surface elem is on - if ( aShapeId != helper.GetSubShapeID() ) { - surface.Nullify(); - TopoDS_Shape shape; - if ( aShapeId > 0 ) - shape = aMesh->IndexToShape( aShapeId ); - if ( !shape.IsNull() && shape.ShapeType() == TopAbs_FACE ) { - TopoDS_Face face = TopoDS::Face( shape ); - surface = BRep_Tool::Surface( face ); - if ( !surface.IsNull() ) - helper.SetSubShape( shape ); - } - } - // find middle point for (0,1,2,3) - // and create a node in this point; - const SMDS_MeshNode* newN = 0; + helper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem )); if ( aNodes.size() == 9 ) { - // SMDSEntity_BiQuad_Quadrangle - newN = aNodes.back(); - } - else - { - gp_XYZ p( 0,0,0 ); - if ( surface.IsNull() ) - { - for ( int i = 0; i < 4; i++ ) - p += gp_XYZ(aNodes[i]->X(), aNodes[i]->Y(), aNodes[i]->Z() ); - p /= 4; - } + helper.SetIsBiQuadratic( true ); + if ( aBadRate1 <= aBadRate2 ) + helper.AddTLinkNode( aNodes[0], aNodes[2], aNodes[8] ); else - { - const SMDS_MeshNode* inFaceNode = 0; - if ( helper.GetNodeUVneedInFaceNode() ) - for ( size_t i = 0; i < aNodes.size() && !inFaceNode; ++i ) - if ( aNodes[ i ]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) - inFaceNode = aNodes[ i ]; - - TopoDS_Face face = TopoDS::Face( helper.GetSubShape() ); - gp_XY uv( 0,0 ); - for ( int i = 0; i < 4; i++ ) - uv += helper.GetNodeUV( face, aNodes[i], inFaceNode ); - uv /= 4.; - p = surface->Value( uv.X(), uv.Y() ).XYZ(); - } - newN = aMesh->AddNode( p.X(), p.Y(), p.Z() ); - myLastCreatedNodes.Append(newN); + helper.AddTLinkNode( aNodes[1], aNodes[3], aNodes[8] ); } // create a new element if ( aBadRate1 <= aBadRate2 ) { - newElem1 = aMesh->AddFace(aNodes[2], aNodes[3], aNodes[0], - aNodes[6], aNodes[7], newN ); - newElem2 = aMesh->AddFace(aNodes[2], aNodes[0], aNodes[1], - newN, aNodes[4], aNodes[5] ); + newElem1 = helper.AddFace( aNodes[2], aNodes[3], aNodes[0] ); + newElem2 = helper.AddFace( aNodes[2], aNodes[0], aNodes[1] ); } else { - newElem1 = aMesh->AddFace(aNodes[3], aNodes[0], aNodes[1], - aNodes[7], aNodes[4], newN ); - newElem2 = aMesh->AddFace(aNodes[3], aNodes[1], aNodes[2], - newN, aNodes[5], aNodes[6] ); + newElem1 = helper.AddFace( aNodes[3], aNodes[0], aNodes[1] ); + newElem2 = helper.AddFace( aNodes[3], aNodes[1], aNodes[2] ); } } // quadratic case @@ -1396,15 +1355,138 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, // put a new triangle on the same shape if ( aShapeId ) - { - aMesh->SetMeshElementOnShape( newElem1, aShapeId ); - aMesh->SetMeshElementOnShape( newElem2, aShapeId ); - } + aMesh->SetMeshElementOnShape( newElem1, aShapeId ); + aMesh->SetMeshElementOnShape( newElem2, aShapeId ); + aMesh->RemoveElement( elem ); } return true; } +//======================================================================= +/*! + * \brief Split each of given quadrangles into 4 triangles. + * \param theElems - The faces to be splitted. If empty all faces are split. + */ +//======================================================================= + +void SMESH_MeshEditor::QuadTo4Tri (TIDSortedElemSet & theElems) +{ + myLastCreatedElems.Clear(); + myLastCreatedNodes.Clear(); + + SMESH_MesherHelper helper( *GetMesh() ); + helper.SetElementsOnShape( true ); + + SMDS_ElemIteratorPtr faceIt; + if ( theElems.empty() ) faceIt = GetMeshDS()->elementsIterator(SMDSAbs_Face); + else faceIt = elemSetIterator( theElems ); + + bool checkUV; + gp_XY uv [9]; uv[8] = gp_XY(0,0); + gp_XYZ xyz[9]; + vector< const SMDS_MeshNode* > nodes; + SMESHDS_SubMesh* subMeshDS; + TopoDS_Face F; + Handle(Geom_Surface) surface; + TopLoc_Location loc; + + while ( faceIt->more() ) + { + const SMDS_MeshElement* quad = faceIt->next(); + if ( !quad || quad->NbCornerNodes() != 4 ) + continue; + + // get a surface the quad is on + + if ( quad->getshapeId() < 1 ) + { + F.Nullify(); + helper.SetSubShape( 0 ); + subMeshDS = 0; + } + else if ( quad->getshapeId() != helper.GetSubShapeID() ) + { + helper.SetSubShape( quad->getshapeId() ); + if ( !helper.GetSubShape().IsNull() && + helper.GetSubShape().ShapeType() == TopAbs_FACE ) + { + F = TopoDS::Face( helper.GetSubShape() ); + surface = BRep_Tool::Surface( F, loc ); + subMeshDS = GetMeshDS()->MeshElements( quad->getshapeId() ); + } + else + { + helper.SetSubShape( 0 ); + subMeshDS = 0; + } + } + + // create a central node + + const SMDS_MeshNode* nCentral; + nodes.assign( quad->begin_nodes(), quad->end_nodes() ); + + if ( nodes.size() == 9 ) + { + nCentral = nodes.back(); + } + else + { + size_t iN = 0; + if ( F.IsNull() ) + { + for ( ; iN < nodes.size(); ++iN ) + xyz[ iN ] = SMESH_TNodeXYZ( nodes[ iN ] ); + + for ( ; iN < 8; ++iN ) // mid-side points of a linear qudrangle + xyz[ iN ] = 0.5 * ( xyz[ iN - 4 ] + xyz[( iN - 3 )%4 ] ); + + xyz[ 8 ] = helper.calcTFI( 0.5, 0.5, + xyz[0], xyz[1], xyz[2], xyz[3], + xyz[4], xyz[5], xyz[6], xyz[7] ); + } + else + { + for ( ; iN < nodes.size(); ++iN ) + uv[ iN ] = helper.GetNodeUV( F, nodes[iN], nodes[(iN+2)%4], &checkUV ); + + for ( ; iN < 8; ++iN ) // UV of mid-side points of a linear qudrangle + uv[ iN ] = helper.GetMiddleUV( surface, uv[ iN - 4 ], uv[( iN - 3 )%4 ] ); + + uv[ 8 ] = helper.calcTFI( 0.5, 0.5, + uv[0], uv[1], uv[2], uv[3], + uv[4], uv[5], uv[6], uv[7] ); + + gp_Pnt p = surface->Value( uv[8].X(), uv[8].Y() ).Transformed( loc ); + xyz[ 8 ] = p.XYZ(); + } + + nCentral = helper.AddNode( xyz[8].X(), xyz[8].Y(), xyz[8].Z(), /*id=*/0, + uv[8].X(), uv[8].Y() ); + myLastCreatedNodes.Append( nCentral ); + } + + // create 4 triangles + + GetMeshDS()->RemoveFreeElement( quad, subMeshDS, /*fromGroups=*/false ); + + helper.SetIsQuadratic ( nodes.size() > 4 ); + helper.SetIsBiQuadratic( nodes.size() == 9 ); + if ( helper.GetIsQuadratic() ) + helper.AddTLinks( static_cast< const SMDS_MeshFace*>( quad )); + + for ( int i = 0; i < 4; ++i ) + { + SMDS_MeshElement* tria = helper.AddFace( nodes[ i ], + nodes[(i+1)%4], + nCentral ); + ReplaceElemInGroups( tria, quad, GetMeshDS() ); + myLastCreatedElems.Append( tria ); + } + } +} + //======================================================================= //function : BestSplit //purpose : Find better diagonal for cutting. @@ -1817,6 +1899,20 @@ namespace }; } // namespace +class TElemToDelete +{ +public: + TElemToDelete(const SMDS_MeshElement* theElem, SMESHDS_SubMesh* theSubMesh) + { + elem = theElem; + subMesh = theSubMesh; + } + const SMDS_MeshElement* Elem() const {return elem;} + SMESHDS_SubMesh* Submesh() {return subMesh;} + const SMDS_MeshElement* elem; + SMESHDS_SubMesh* subMesh; +}; + //======================================================================= //function : SplitVolumesIntoTetra //purpose : Split volume elements into tetrahedra. @@ -1842,6 +1938,7 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, double bc[3]; TIDSortedElemSet::const_iterator elem = theElems.begin(); + std::vector elem_to_delete; for ( ; elem != theElems.end(); ++elem ) { if ( (*elem)->GetType() != SMDSAbs_Volume ) @@ -2005,11 +2102,26 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } ReplaceElemInGroups( face, triangles, GetMeshDS() ); GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false ); +// TElemToDelete faceToDelete(face, fSubMesh); +// elem_to_delete.push_back(faceToDelete); } } // loop on volume faces to split them into triangles - GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false ); +// GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false ); + + // rnc : don't delete the elem here because it results in a mesh with a free + // ID at the beginning of the ID list. The first tetra is then inserted in O(1) + // but the second one is inserted in O(n), then the whole procedure has almost a O(n^2) + // complexity. If all elements to remove are stored and removed after tetra creation + // we get a O(n) complexity for the whole procedure. + // The memory cost is at worst a 6*n*constant memory occupation (where n is the number of elements) + // before deletion of the hexas and then 5*n*constant instead of a maximum of 5*n*constant. + // So there is a transient 1/5*(memory occupation) additional cost. + + // Store the elements to delete + TElemToDelete elemToDelete(*elem, subMesh); + elem_to_delete.push_back(elemToDelete); if ( geomType == SMDSEntity_TriQuad_Hexa ) { @@ -2020,6 +2132,13 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } } // loop on volumes to split + // Delete stored elements + std::vector::iterator it; + for( it = elem_to_delete.begin(); it!= elem_to_delete.end(); it++) + { + GetMeshDS()->RemoveFreeElement( it->Elem(), it->Submesh(), /*fromGroups=*/false ); + } + myLastCreatedNodes = newNodes; myLastCreatedElems = newElems; } @@ -2380,18 +2499,19 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems, map< const SMDS_MeshElement*, set< SMESH_TLink > >::iterator itEL; TIDSortedElemSet::iterator itElem; - for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) { + for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) + { const SMDS_MeshElement* elem = *itElem; if(!elem || elem->GetType() != SMDSAbs_Face ) continue; - bool IsTria = elem->NbNodes()==3 || (elem->NbNodes()==6 && elem->IsQuadratic()); - if(!IsTria) continue; + bool IsTria = ( elem->NbCornerNodes()==3 ); + if (!IsTria) continue; // retrieve element nodes const SMDS_MeshNode* aNodes [4]; - SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + SMDS_NodeIteratorPtr itN = elem->nodeIterator(); int i = 0; - while ( i<3 ) - aNodes[ i++ ] = cast2Node( itN->next() ); + while ( i < 3 ) + aNodes[ i++ ] = itN->next(); aNodes[ 3 ] = aNodes[ 0 ]; // fill maps @@ -2536,31 +2656,31 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems, } // Make quadrangles - // and remove fused elems and removed links from the maps + // and remove fused elems and remove links from the maps mapEl_setLi.erase( tr1 ); - if ( Ok12 ) { + if ( Ok12 ) + { mapEl_setLi.erase( tr2 ); mapLi_listEl.erase( *link12 ); - if(tr1->NbNodes()==3) { + if ( tr1->NbNodes() == 3 ) + { const SMDS_MeshElement* newElem = 0; newElem = aMesh->AddFace(n12[0], n12[1], n12[2], n12[3] ); myLastCreatedElems.Append(newElem); AddToSameGroups( newElem, tr1, aMesh ); int aShapeId = tr1->getshapeId(); if ( aShapeId ) - { - aMesh->SetMeshElementOnShape( newElem, aShapeId ); - } + aMesh->SetMeshElementOnShape( newElem, aShapeId ); aMesh->RemoveElement( tr1 ); aMesh->RemoveElement( tr2 ); } else { - const SMDS_MeshNode* N1 [6]; - const SMDS_MeshNode* N2 [6]; - GetNodesFromTwoTria(tr1,tr2,N1,N2); - // now we receive following N1 and N2 (using numeration as above image) + vector< const SMDS_MeshNode* > N1; + vector< const SMDS_MeshNode* > N2; + getNodesFromTwoTria(tr1,tr2,N1,N2); + // now we receive following N1 and N2 (using numeration as in image in InverseDiag()) // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6) - // i.e. first nodes from both arrays determ new diagonal + // i.e. first nodes from both arrays form a new diagonal const SMDS_MeshNode* aNodes[8]; aNodes[0] = N1[0]; aNodes[1] = N1[1]; @@ -2571,44 +2691,50 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems, aNodes[6] = N2[3]; aNodes[7] = N1[5]; const SMDS_MeshElement* newElem = 0; - newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], - aNodes[4], aNodes[5], aNodes[6], aNodes[7]); + if ( N1.size() == 7 || N2.size() == 7 ) // biquadratic + newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], + aNodes[4], aNodes[5], aNodes[6], aNodes[7], N1[4]); + else + newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], + aNodes[4], aNodes[5], aNodes[6], aNodes[7]); myLastCreatedElems.Append(newElem); AddToSameGroups( newElem, tr1, aMesh ); int aShapeId = tr1->getshapeId(); if ( aShapeId ) - { - aMesh->SetMeshElementOnShape( newElem, aShapeId ); - } + aMesh->SetMeshElementOnShape( newElem, aShapeId ); aMesh->RemoveElement( tr1 ); aMesh->RemoveElement( tr2 ); // remove middle node (9) - GetMeshDS()->RemoveNode( N1[4] ); + if ( N1[4]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N1[4] ); + if ( N1.size() == 7 && N1[6]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N1[6] ); + if ( N2.size() == 7 && N2[6]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N2[6] ); } } - else if ( Ok13 ) { + else if ( Ok13 ) + { mapEl_setLi.erase( tr3 ); mapLi_listEl.erase( *link13 ); - if(tr1->NbNodes()==3) { + if ( tr1->NbNodes() == 3 ) { const SMDS_MeshElement* newElem = 0; newElem = aMesh->AddFace(n13[0], n13[1], n13[2], n13[3] ); myLastCreatedElems.Append(newElem); AddToSameGroups( newElem, tr1, aMesh ); int aShapeId = tr1->getshapeId(); if ( aShapeId ) - { - aMesh->SetMeshElementOnShape( newElem, aShapeId ); - } + aMesh->SetMeshElementOnShape( newElem, aShapeId ); aMesh->RemoveElement( tr1 ); aMesh->RemoveElement( tr3 ); } else { - const SMDS_MeshNode* N1 [6]; - const SMDS_MeshNode* N2 [6]; - GetNodesFromTwoTria(tr1,tr3,N1,N2); + vector< const SMDS_MeshNode* > N1; + vector< const SMDS_MeshNode* > N2; + getNodesFromTwoTria(tr1,tr3,N1,N2); // now we receive following N1 and N2 (using numeration as above image) // tria1 : (1 2 4 5 9 7) and tria2 : (3 4 2 8 9 6) - // i.e. first nodes from both arrays determ new diagonal + // i.e. first nodes from both arrays form a new diagonal const SMDS_MeshNode* aNodes[8]; aNodes[0] = N1[0]; aNodes[1] = N1[1]; @@ -2619,19 +2745,26 @@ bool SMESH_MeshEditor::TriToQuad (TIDSortedElemSet & theElems, aNodes[6] = N2[3]; aNodes[7] = N1[5]; const SMDS_MeshElement* newElem = 0; - newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], - aNodes[4], aNodes[5], aNodes[6], aNodes[7]); + if ( N1.size() == 7 || N2.size() == 7 ) // biquadratic + newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], + aNodes[4], aNodes[5], aNodes[6], aNodes[7], N1[4]); + else + newElem = aMesh->AddFace(aNodes[0], aNodes[1], aNodes[2], aNodes[3], + aNodes[4], aNodes[5], aNodes[6], aNodes[7]); myLastCreatedElems.Append(newElem); AddToSameGroups( newElem, tr1, aMesh ); int aShapeId = tr1->getshapeId(); if ( aShapeId ) - { - aMesh->SetMeshElementOnShape( newElem, aShapeId ); - } + aMesh->SetMeshElementOnShape( newElem, aShapeId ); aMesh->RemoveElement( tr1 ); aMesh->RemoveElement( tr3 ); // remove middle node (9) - GetMeshDS()->RemoveNode( N1[4] ); + if ( N1[4]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N1[4] ); + if ( N1.size() == 7 && N1[6]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N1[6] ); + if ( N2.size() == 7 && N2[6]->NbInverseElements() == 0 ) + aMesh->RemoveNode( N2[6] ); } } @@ -3590,39 +3723,34 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet & theElems, if ( isQuadratic ) { SMESH_MesherHelper helper( *GetMesh() ); - if ( !face.IsNull() ) - helper.SetSubShape( face ); + helper.SetSubShape( face ); + vector nodes; + bool checkUV; list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin(); - for ( ; elemIt != elemsOnFace.end(); ++elemIt ) { - const SMDS_VtkFace* QF = - dynamic_cast (*elemIt); - if(QF && QF->IsQuadratic()) { - vector Ns; - Ns.reserve(QF->NbNodes()+1); - SMDS_ElemIteratorPtr anIter = QF->interlacedNodesElemIterator(); - while ( anIter->more() ) - Ns.push_back( cast2Node(anIter->next()) ); - Ns.push_back( Ns[0] ); - double x, y, z; - for(int i=0; iNbNodes(); i=i+2) { - if ( !surface.IsNull() ) { - gp_XY uv1 = helper.GetNodeUV( face, Ns[i], Ns[i+2] ); - gp_XY uv2 = helper.GetNodeUV( face, Ns[i+2], Ns[i] ); - gp_XY uv = ( uv1 + uv2 ) / 2.; - gp_Pnt xyz = surface->Value( uv.X(), uv.Y() ); - x = xyz.X(); y = xyz.Y(); z = xyz.Z(); + for ( ; elemIt != elemsOnFace.end(); ++elemIt ) + { + const SMDS_MeshElement* QF = *elemIt; + if ( QF->IsQuadratic() ) + { + nodes.assign( SMDS_MeshElement::iterator( QF->interlacedNodesElemIterator() ), + SMDS_MeshElement::iterator() ); + nodes.push_back( nodes[0] ); + gp_Pnt xyz; + for (size_t i = 1; i < nodes.size(); i += 2 ) // i points to a medium node + { + if ( !surface.IsNull() ) + { + gp_XY uv1 = helper.GetNodeUV( face, nodes[i-1], nodes[i+1], &checkUV ); + gp_XY uv2 = helper.GetNodeUV( face, nodes[i+1], nodes[i-1], &checkUV ); + gp_XY uv = helper.GetMiddleUV( surface, uv1, uv2 ); + xyz = surface->Value( uv.X(), uv.Y() ); } else { - x = (Ns[i]->X() + Ns[i+2]->X())/2; - y = (Ns[i]->Y() + Ns[i+2]->Y())/2; - z = (Ns[i]->Z() + Ns[i+2]->Z())/2; - } - if( fabs( Ns[i+1]->X() - x ) > disttol || - fabs( Ns[i+1]->Y() - y ) > disttol || - fabs( Ns[i+1]->Z() - z ) > disttol ) { - // we have to move i+1 node - aMesh->MoveNode( Ns[i+1], x, y, z ); + xyz = 0.5 * ( SMESH_TNodeXYZ( nodes[i-1] ) + SMESH_TNodeXYZ( nodes[i+1] )); } + if (( SMESH_TNodeXYZ( nodes[i] ) - xyz.XYZ() ).Modulus() > disttol ) + // we have to move a medium node + aMesh->MoveNode( nodes[i], xyz.X(), xyz.Y(), xyz.Z() ); } } } @@ -3649,7 +3777,7 @@ static bool isReverse(const SMDS_MeshElement* face, SMESH_TNodeXYZ pP = prevNodes[ iNotSame ]; SMESH_TNodeXYZ pN = nextNodes[ iNotSame ]; gp_XYZ extrDir( pN - pP ), faceNorm; - SMESH_Algo::FaceNormal( face, faceNorm, /*normalized=*/false ); + SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false ); return faceNorm * extrDir < 0.0; } @@ -3904,7 +4032,8 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem, } break; } - case SMDSEntity_Quad_Triangle: { // sweep Quadratic TRIANGLE ---> + case SMDSEntity_Quad_Triangle: // sweep (Bi)Quadratic TRIANGLE ---> + case SMDSEntity_BiQuad_Triangle: /* ??? */ { if ( nbDouble+nbSame != 3 ) break; if(nbSame==0) { // ---> pentahedron with 15 nodes @@ -4117,7 +4246,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, ASSERT( newElemsMap.size() == elemNewNodesMap.size() ); SMESHDS_Mesh* aMesh = GetMeshDS(); - // Find nodes belonging to only one initial element - sweep them to get edges. + // Find nodes belonging to only one initial element - sweep them into edges. TNodeOfNodeListMapItr nList = mapNewNodes.begin(); for ( ; nList != mapNewNodes.end(); nList++ ) @@ -4204,7 +4333,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, const SMDS_MeshNode* n1 = vecNewNodes[ iNode ]->first; const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first; // check if a link n1-n2 is free - if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet )) { + if ( ! SMESH_MeshAlgos::FindFaceInSet ( n1, n2, elemSet, avoidSet )) { hasFreeLinks = true; // make a new edge and a ceiling for a new edge const SMDS_MeshElement* edge; @@ -4230,9 +4359,9 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, const SMDS_MeshNode* n2 = vecNewNodes[ iNext ]->first; const SMDS_MeshNode* n3 = vecNewNodes[ iNode+nbn ]->first; // check if a link is free - if ( ! SMESH_MeshEditor::FindFaceInSet ( n1, n2, elemSet, avoidSet ) && - ! SMESH_MeshEditor::FindFaceInSet ( n1, n3, elemSet, avoidSet ) && - ! SMESH_MeshEditor::FindFaceInSet ( n3, n2, elemSet, avoidSet ) ) { + if ( ! SMESH_MeshAlgos::FindFaceInSet ( n1, n2, elemSet, avoidSet ) && + ! SMESH_MeshAlgos::FindFaceInSet ( n1, n3, elemSet, avoidSet ) && + ! SMESH_MeshAlgos::FindFaceInSet ( n3, n2, elemSet, avoidSet ) ) { hasFreeLinks = true; // make an edge and a ceiling for a new edge // find medium node @@ -4261,10 +4390,15 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, int iVol, volNb, nbVolumesByStep = newVolumes.size() / nbSteps; set initNodeSet, topNodeSet, faceNodeSet; + set initNodeSetNoCenter/*, topNodeSetNoCenter*/; for ( iNode = 0; iNode < nbNodes; iNode++ ) { initNodeSet.insert( vecNewNodes[ iNode ]->first ); topNodeSet .insert( vecNewNodes[ iNode ]->second.back() ); } + if ( isQuadratic && nbNodes % 2 ) { // node set for the case of a biquadratic + initNodeSetNoCenter = initNodeSet; // swept face and a not biquadratic volume + initNodeSetNoCenter.erase( vecNewNodes.back()->first ); + } for ( volNb = 0; volNb < nbVolumesByStep; volNb++ ) { list::iterator v = newVolumes.begin(); std::advance( v, volNb ); @@ -4280,6 +4414,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, { if ( nbSteps == 1 && faceNodeSet == topNodeSet ) continue; + if ( faceNodeSet == initNodeSetNoCenter ) + continue; freeInd.push_back( iF ); // find source edge of a free face iF vector commonNodes; // shared by the initial and free faces @@ -4452,56 +4588,48 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, // Make a ceiling face with a normal external to a volume + // use SMDS_VolumeTool to get a correctly ordered nodes of a ceiling face SMDS_VolumeTool lastVol( itElem->second.back(), /*ignoreCentralNodes=*/false ); - int iF = lastVol.GetFaceIndex( aFaceLastNodes ); + + if ( iF < 0 && isQuadratic && nbNodes % 2 ) { // remove a central node of biquadratic + aFaceLastNodes.erase( vecNewNodes.back()->second.back() ); + iF = lastVol.GetFaceIndex( aFaceLastNodes ); + } if ( iF >= 0 ) { lastVol.SetExternalNormal(); const SMDS_MeshNode** nodes = lastVol.GetFaceNodes( iF ); int nbn = lastVol.NbFaceNodes( iF ); - if ( nbn == 3 ) { - if (!hasFreeLinks || - !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ])) - myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ] )); - } - else if ( nbn == 4 ) - { - if (!hasFreeLinks || - !aMesh->FindFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ])) - myLastCreatedElems.Append(aMesh->AddFace( nodes[ 0 ], nodes[ 1 ], nodes[ 2 ], nodes[ 3 ])); - } - else if ( nbn == 6 && isQuadratic ) - { - if (!hasFreeLinks || - !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[1], nodes[3], nodes[5]) ) - myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], - nodes[1], nodes[3], nodes[5])); - } - else if ( nbn == 8 && isQuadratic ) + // we do not use this->AddElement() because nodes are interlaced + vector nodeVec( nodes, nodes+nbn ); + if ( !hasFreeLinks || + !aMesh->FindElement( nodeVec, SMDSAbs_Face, /*noMedium=*/false) ) { - if (!hasFreeLinks || - !aMesh->FindFace(nodes[0], nodes[2], nodes[4], nodes[6], - nodes[1], nodes[3], nodes[5], nodes[7]) ) - myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6], - nodes[1], nodes[3], nodes[5], nodes[7])); - } - else if ( nbn == 9 && isQuadratic ) - { - if (!hasFreeLinks || - !aMesh->FindElement(vector( nodes, nodes+nbn ), - SMDSAbs_Face, /*noMedium=*/false) ) - myLastCreatedElems.Append(aMesh->AddFace(nodes[0], nodes[2], nodes[4], nodes[6], - nodes[1], nodes[3], nodes[5], nodes[7], - nodes[8])); - } - else { - vector polygon_nodes ( nodes, nodes + nbn ); - if (!hasFreeLinks || !aMesh->FindFace(polygon_nodes)) - myLastCreatedElems.Append(aMesh->AddPolygonalFace(polygon_nodes)); - } + if ( nbn == 3 ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[1], nodes[2] )); + + else if ( nbn == 4 ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[1], nodes[2], nodes[3])); + + else if ( nbn == 6 && isQuadratic ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[2], nodes[4], + nodes[1], nodes[3], nodes[5])); + else if ( nbn == 7 && isQuadratic ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[2], nodes[4], + nodes[1], nodes[3], nodes[5], nodes[6])); + else if ( nbn == 8 && isQuadratic ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[2], nodes[4], nodes[6], + nodes[1], nodes[3], nodes[5], nodes[7])); + else if ( nbn == 9 && isQuadratic ) + myLastCreatedElems.Append(aMesh->AddFace( nodes[0], nodes[2], nodes[4], nodes[6], + nodes[1], nodes[3], nodes[5], nodes[7], + nodes[8])); + else + myLastCreatedElems.Append(aMesh->AddPolygonalFace( nodeVec )); - while ( srcElements.Length() < myLastCreatedElems.Length() ) - srcElements.Append( elem ); + while ( srcElements.Length() < myLastCreatedElems.Length() ) + srcElements.Append( elem ); + } } } // loop on swept elements } @@ -5053,7 +5181,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, const TopoDS_Shape& aS = theTrack->GetShapeToMesh(); - if( aS == SMESH_Mesh::PseudoShape() ) { + if ( !theTrack->HasShapeToMesh() ) { //Mesh without shape const SMDS_MeshNode* currentNode = NULL; const SMDS_MeshNode* prevNode = theN1; @@ -5176,16 +5304,8 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, if ( BRep_Tool::Degenerated( aTrackEdge ) ) return EXTR_BAD_PATH_SHAPE; TopExp::Vertices( aTrackEdge, aV1, aV2 ); - const SMDS_MeshNode* aN1 = 0; - const SMDS_MeshNode* aN2 = 0; - if ( theTrack->GetSubMesh( aV1 ) && theTrack->GetSubMesh( aV1 )->GetSubMeshDS() ) { - aItN = theTrack->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes(); - aN1 = aItN->next(); - } - if ( theTrack->GetSubMesh( aV2 ) && theTrack->GetSubMesh( aV2 )->GetSubMeshDS() ) { - aItN = theTrack->GetSubMesh( aV2 )->GetSubMeshDS()->GetNodes(); - aN2 = aItN->next(); - } + const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS ); + const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS ); // starting node must be aN1 or aN2 if ( !( aN1 == theN1 || aN2 == theN1 ) ) return EXTR_BAD_STARTING_NODE; @@ -5233,17 +5353,8 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, if ( aVprev.IsNull() ) { // if previous vertex is not yet defined, it means that we in the beginning of wire // and we have to find initial vertex corresponding to starting node theN1 - const SMDS_MeshNode* aN1 = 0; - const SMDS_MeshNode* aN2 = 0; - - if ( locTrack->GetFather()->GetSubMesh(aV1) && locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS() ) { - aItN = locTrack->GetFather()->GetSubMesh(aV1)->GetSubMeshDS()->GetNodes(); - aN1 = aItN->next(); - } - if ( locTrack->GetFather()->GetSubMesh(aV2) && locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS() ) { - aItN = locTrack->GetFather()->GetSubMesh(aV2)->GetSubMeshDS()->GetNodes(); - aN2 = aItN->next(); - } + const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS ); + const SMDS_MeshNode* aN2 = SMESH_Algo::VertexNode( aV2, pMeshDS ); // starting node must be aN1 or aN2 aN1isOK = ( aN1 && aN1 == theN1 ); aN2isOK = ( aN2 && aN2 == theN1 ); @@ -5277,27 +5388,21 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, } } list< list >::iterator itLLPP = LLPPs.begin(); - list firstList = *itLLPP; - list::iterator itPP = firstList.begin(); - for(; itPP!=firstList.end(); itPP++) { - fullList.push_back( *itPP ); - } + list& firstList = *itLLPP; + fullList.splice( fullList.end(), firstList ); + SMESH_MeshEditor_PathPoint PP1 = fullList.back(); fullList.pop_back(); itLLPP++; for(; itLLPP!=LLPPs.end(); itLLPP++) { - list currList = *itLLPP; - itPP = currList.begin(); + list& currList = *itLLPP; SMESH_MeshEditor_PathPoint PP2 = currList.front(); gp_Dir D1 = PP1.Tangent(); gp_Dir D2 = PP2.Tangent(); gp_Dir Dnew( ( D1.XYZ() + D2.XYZ() ) / 2 ); PP1.SetTangent(Dnew); fullList.push_back(PP1); - itPP++; - for(; itPP!=currList.end(); itPP++) { - fullList.push_back( *itPP ); - } + fullList.splice( fullList.end(), currList, ++currList.begin(), currList.end() ); PP1 = fullList.back(); fullList.pop_back(); } @@ -5319,9 +5424,9 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, //purpose : auxilary for ExtrusionAlongTrack //======================================================================= SMESH_MeshEditor::Extrusion_Error -SMESH_MeshEditor::MakeEdgePathPoints(std::list& aPrms, - const TopoDS_Edge& aTrackEdge, - bool FirstIsStart, +SMESH_MeshEditor::MakeEdgePathPoints(std::list& aPrms, + const TopoDS_Edge& aTrackEdge, + bool FirstIsStart, list& LPP) { Standard_Real aTx1, aTx2, aL2, aTolVec, aTolVec2; @@ -5374,63 +5479,43 @@ SMESH_MeshEditor::MakeEdgePathPoints(std::list& aPrms, //purpose : auxilary for ExtrusionAlongTrack //======================================================================= SMESH_MeshEditor::Extrusion_Error -SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, +SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, list& fullList, - const bool theHasAngles, - list& theAngles, - const bool theLinearVariation, - const bool theHasRefPoint, - const gp_Pnt& theRefPoint, - const bool theMakeGroups) + const bool theHasAngles, + list& theAngles, + const bool theLinearVariation, + const bool theHasRefPoint, + const gp_Pnt& theRefPoint, + const bool theMakeGroups) { - MESSAGE("MakeExtrElements"); - //cout<<"MakeExtrElements fullList.size() = "< aPPs(aNbTP); + const int aNbTP = fullList.size(); // Angles - if( theHasAngles && theAngles.size()>0 && theLinearVariation ) { + if( theHasAngles && !theAngles.empty() && theLinearVariation ) LinearAngleVariation(aNbTP-1, theAngles); - } - vector aAngles( aNbTP ); - int j = 0; - for(; j::iterator aItD = theAngles.begin(); - for ( j=1; (aItD != theAngles.end()) && (j aPPs; list::iterator itPP = fullList.begin(); - for(; itPP!=fullList.end(); itPP++) { - j++; - SMESH_MeshEditor_PathPoint PP = *itPP; - PP.SetAngle(aAngles[j]); - aPPs[j] = PP; + list::iterator itAngles = theAngles.begin(); + aPPs.push_back( *itPP++ ); + for( ; itPP != fullList.end(); itPP++) { + aPPs.push_back( *itPP ); + if ( theHasAngles && itAngles != theAngles.end() ) + aPPs.back().SetAngle( *itAngles ); } - TNodeOfNodeListMap mapNewNodes; + TNodeOfNodeListMap mapNewNodes; TElemOfVecOfNnlmiMap mapElemNewNodes; - TElemOfElemListMap newElemsMap; + TElemOfElemListMap newElemsMap; TIDSortedElemSet::iterator itElem; - double aX, aY, aZ; - int aNb; - SMDSAbs_ElementType aTypeE; // source elements for each generated one SMESH_SequenceOfElemPtr srcElems, srcNodes; // 3. Center of rotation aV0 gp_Pnt aV0 = theRefPoint; - gp_XYZ aGC; - if ( !theHasRefPoint ) { - aNb = 0; - aGC.SetCoord( 0.,0.,0. ); + if ( !theHasRefPoint ) + { + gp_XYZ aGC( 0.,0.,0. ); + TIDSortedElemSet newNodes; itElem = theElements.begin(); for ( ; itElem != theElements.end(); itElem++ ) { @@ -5438,25 +5523,14 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, SMDS_ElemIteratorPtr itN = elem->nodesIterator(); while ( itN->more() ) { - const SMDS_MeshNode* node = static_cast( itN->next() ); - aX = node->X(); - aY = node->Y(); - aZ = node->Z(); - - if ( mapNewNodes.find( node ) == mapNewNodes.end() ) { - list aLNx; - mapNewNodes[node] = aLNx; - // - gp_XYZ aXYZ( aX, aY, aZ ); - aGC += aXYZ; - ++aNb; - } + const SMDS_MeshElement* node = itN->next(); + if ( newNodes.insert( node ).second ) + aGC += SMESH_TNodeXYZ( node ); } } - aGC /= aNb; + aGC /= newNodes.size(); aV0.SetXYZ( aGC ); } // if (!theHasRefPoint) { - mapNewNodes.clear(); // 4. Processing the elements SMESHDS_Mesh* aMesh = GetMeshDS(); @@ -5464,7 +5538,7 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, for ( itElem = theElements.begin(); itElem != theElements.end(); itElem++ ) { // check element type const SMDS_MeshElement* elem = *itElem; - aTypeE = elem->GetType(); + SMDSAbs_ElementType aTypeE = elem->GetType(); if ( !elem || ( aTypeE != SMDSAbs_Face && aTypeE != SMDSAbs_Edge ) ) continue; @@ -5486,8 +5560,6 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, list& listNewNodes = nIt->second; // make new nodes - aX = node->X(); aY = node->Y(); aZ = node->Z(); - Standard_Real aAngle1x, aAngleT1T0, aTolAng; gp_Pnt aP0x, aP1x, aPN0, aPN1, aV0x, aV1x; gp_Ax1 anAx1, anAxT1T0; @@ -5496,17 +5568,17 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements, aTolAng=1.e-4; aV0x = aV0; - aPN0.SetCoord(aX, aY, aZ); + aPN0 = SMESH_TNodeXYZ( node ); const SMESH_MeshEditor_PathPoint& aPP0 = aPPs[0]; aP0x = aPP0.Pnt(); aDT0x= aPP0.Tangent(); //cout<<"j = 0 PP: Pnt("<AddNode( aX, aY, aZ ); + const SMDS_MeshNode* newNode = aMesh->AddNode( aPN1.X(), aPN1.Y(), aPN1.Z() ); myLastCreatedNodes.Append(newNode); srcNodes.Append( node ); listNewNodes.push_back( newNode ); @@ -6028,8 +6097,8 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, const SMDSAbs_ElementType type = groupDS->GetType(); SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type ); SMESHDS_Group* newTopGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), type ); - groupsByType[ groupDS->GetType() ].push_back( make_tuple( groupDS, newGroup, newTopGroup )); - orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() ); + groupsByType[ type ].push_back( make_tuple( groupDS, newGroup, newTopGroup )); + orderedOldNewGroups.push_back( & groupsByType[ type ].back() ); } // Loop on nodes and elements to add them in new groups @@ -6065,1488 +6134,136 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, if ( resElem != sourceElem ) resultElems.push_back( resElem ); - // add resultElems to groups made by ones the sourceElem belongs to - list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end(); - for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew ) - { - SMESHDS_GroupBase* oldGroup = gOldNew->get<0>(); - if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup - { - // fill in a new group - SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup(); - list< const SMDS_MeshElement* > rejectedElems; // elements of other type - list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt; - for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt ) - if ( !newGroup.Add( *resElemIt )) - rejectedElems.push_back( *resElemIt ); - - // fill "top" group - if ( !rejectedElems.empty() ) - { - SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup(); - resLast = rejectedElems.end(); - for ( resElemIt = rejectedElems.begin(); resElemIt != resLast; ++resElemIt ) - !newTopGroup.Add( *resElemIt ); - } - } - } - } // loop on created elements - }// loop on nodes and elements - - // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups - - list topGrouIds; - for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i ) - { - SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->get<0>(); - SMESHDS_Group* newGroups[2] = { orderedOldNewGroups[i]->get<1>(), - orderedOldNewGroups[i]->get<2>() }; - const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty(); - for ( int is2nd = 0; is2nd < 2; ++is2nd ) - { - SMESHDS_Group* newGroupDS = newGroups[ is2nd ]; - if ( newGroupDS->IsEmpty() ) - { - mesh->GetMeshDS()->RemoveGroup( newGroupDS ); - } - else - { - // set group type - newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() ); - - // make a name - const bool isTop = ( nbNewGroups == 2 && - newGroupDS->GetType() == oldGroupDS->GetType() ); - string name = oldGroupDS->GetStoreName(); - if ( !targetMesh ) { - string suffix = ( isTop ? "top": postfix.c_str() ); - name += "_"; - name += suffix; - int nb = 1; - while ( !groupNames.insert( name ).second ) // name exists - name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << suffix << "_" << nb++; - } - else if ( isTop ) { - name += "_top"; - } - newGroupDS->SetStoreName( name.c_str() ); - - // make a SMESH_Groups - mesh->AddGroup( newGroupDS ); - if ( isTop ) - topGrouIds.push_back( newGroupDS->GetID() ); - else - newGroupIDs->push_back( newGroupDS->GetID() ); - } - } - } - newGroupIDs->splice( newGroupIDs->end(), topGrouIds ); - - return newGroupIDs; -} - -//================================================================================ -/*! - * \brief Return list of group of nodes close to each other within theTolerance - * Search among theNodes or in the whole mesh if theNodes is empty using - * an Octree algorithm - */ -//================================================================================ - -void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes, - const double theTolerance, - TListOfListOfNodes & theGroupsOfNodes) -{ - myLastCreatedElems.Clear(); - myLastCreatedNodes.Clear(); - - if ( theNodes.empty() ) - { // get all nodes in the mesh - SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true); - while ( nIt->more() ) - theNodes.insert( theNodes.end(),nIt->next()); - } - - SMESH_OctreeNode::FindCoincidentNodes ( theNodes, &theGroupsOfNodes, theTolerance); -} - - -//======================================================================= -/*! - * \brief Implementation of search for the node closest to point - */ -//======================================================================= - -struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher -{ - //--------------------------------------------------------------------- - /*! - * \brief Constructor - */ - SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh ) - { - myMesh = ( SMESHDS_Mesh* ) theMesh; - - TIDSortedNodeSet nodes; - if ( theMesh ) { - SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true); - while ( nIt->more() ) - nodes.insert( nodes.end(), nIt->next() ); - } - myOctreeNode = new SMESH_OctreeNode(nodes) ; - - // get max size of a leaf box - SMESH_OctreeNode* tree = myOctreeNode; - while ( !tree->isLeaf() ) - { - SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator(); - if ( cIt->more() ) - tree = cIt->next(); - } - myHalfLeafSize = tree->maxSize() / 2.; - } - - //--------------------------------------------------------------------- - /*! - * \brief Move node and update myOctreeNode accordingly - */ - void MoveNode( const SMDS_MeshNode* node, const gp_Pnt& toPnt ) - { - myOctreeNode->UpdateByMoveNode( node, toPnt ); - myMesh->MoveNode( node, toPnt.X(), toPnt.Y(), toPnt.Z() ); - } - - //--------------------------------------------------------------------- - /*! - * \brief Do it's job - */ - const SMDS_MeshNode* FindClosestTo( const gp_Pnt& thePnt ) - { - map dist2Nodes; - myOctreeNode->NodesAround( thePnt.Coord(), dist2Nodes, myHalfLeafSize ); - if ( !dist2Nodes.empty() ) - return dist2Nodes.begin()->second; - list nodes; - //myOctreeNode->NodesAround( &tgtNode, &nodes, myHalfLeafSize ); - - double minSqDist = DBL_MAX; - if ( nodes.empty() ) // get all nodes of OctreeNode's closest to thePnt - { - // sort leafs by their distance from thePnt - typedef map< double, SMESH_OctreeNode* > TDistTreeMap; - TDistTreeMap treeMap; - list< SMESH_OctreeNode* > treeList; - list< SMESH_OctreeNode* >::iterator trIt; - treeList.push_back( myOctreeNode ); - - gp_XYZ pointNode( thePnt.X(), thePnt.Y(), thePnt.Z() ); - bool pointInside = myOctreeNode->isInside( pointNode, myHalfLeafSize ); - for ( trIt = treeList.begin(); trIt != treeList.end(); ++trIt) - { - SMESH_OctreeNode* tree = *trIt; - if ( !tree->isLeaf() ) // put children to the queue - { - if ( pointInside && !tree->isInside( pointNode, myHalfLeafSize )) continue; - SMESH_OctreeNodeIteratorPtr cIt = tree->GetChildrenIterator(); - while ( cIt->more() ) - treeList.push_back( cIt->next() ); - } - else if ( tree->NbNodes() ) // put a tree to the treeMap - { - const Bnd_B3d& box = *tree->getBox(); - double sqDist = thePnt.SquareDistance( 0.5 * ( box.CornerMin() + box.CornerMax() )); - pair it_in = treeMap.insert( make_pair( sqDist, tree )); - if ( !it_in.second ) // not unique distance to box center - treeMap.insert( it_in.first, make_pair( sqDist + 1e-13*treeMap.size(), tree )); - } - } - // find distance after which there is no sense to check tree's - double sqLimit = DBL_MAX; - TDistTreeMap::iterator sqDist_tree = treeMap.begin(); - if ( treeMap.size() > 5 ) { - SMESH_OctreeNode* closestTree = sqDist_tree->second; - const Bnd_B3d& box = *closestTree->getBox(); - double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() ); - sqLimit = limit * limit; - } - // get all nodes from trees - for ( ; sqDist_tree != treeMap.end(); ++sqDist_tree) { - if ( sqDist_tree->first > sqLimit ) - break; - SMESH_OctreeNode* tree = sqDist_tree->second; - tree->NodesAround( tree->GetNodeIterator()->next(), &nodes ); - } - } - // find closest among nodes - minSqDist = DBL_MAX; - const SMDS_MeshNode* closestNode = 0; - list::iterator nIt = nodes.begin(); - for ( ; nIt != nodes.end(); ++nIt ) { - double sqDist = thePnt.SquareDistance( SMESH_TNodeXYZ( *nIt ) ); - if ( minSqDist > sqDist ) { - closestNode = *nIt; - minSqDist = sqDist; - } - } - return closestNode; - } - - //--------------------------------------------------------------------- - /*! - * \brief Destructor - */ - ~SMESH_NodeSearcherImpl() { delete myOctreeNode; } - - //--------------------------------------------------------------------- - /*! - * \brief Return the node tree - */ - const SMESH_OctreeNode* getTree() const { return myOctreeNode; } - -private: - SMESH_OctreeNode* myOctreeNode; - SMESHDS_Mesh* myMesh; - double myHalfLeafSize; // max size of a leaf box -}; - -//======================================================================= -/*! - * \brief Return SMESH_NodeSearcher - */ -//======================================================================= - -SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher() -{ - return new SMESH_NodeSearcherImpl( GetMeshDS() ); -} - -// ======================================================================== -namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() -{ - const int MaxNbElemsInLeaf = 10; // maximal number of elements in a leaf of tree - const int MaxLevel = 7; // maximal tree height -> nb terminal boxes: 8^7 = 2097152 - const double NodeRadius = 1e-9; // to enlarge bnd box of element - - //======================================================================= - /*! - * \brief Octal tree of bounding boxes of elements - */ - //======================================================================= - - class ElementBndBoxTree : public SMESH_Octree - { - 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 ); - 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():_size(0) {} - SMESH_Octree* newChild() const { return new ElementBndBoxTree; } - void buildChildrenData(); - Bnd_B3d* buildRootBox(); - private: - //!< Bounding box of element - struct ElementBox : public Bnd_B3d - { - const SMDS_MeshElement* _element; - int _refCount; // an ElementBox can be included in several tree branches - ElementBox(const SMDS_MeshElement* elem, double tolerance); - }; - vector< ElementBox* > _elements; - size_t _size; - }; - - //================================================================================ - /*! - * \brief ElementBndBoxTree creation - */ - //================================================================================ - - ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance) - :SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. )) - { - int nbElems = mesh.GetMeshInfo().NbElements( elemType ); - _elements.reserve( nbElems ); - - SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); - while ( elemIt->more() ) - _elements.push_back( new ElementBox( elemIt->next(),tolerance )); - - compute(); - } - - //================================================================================ - /*! - * \brief Destructor - */ - //================================================================================ - - ElementBndBoxTree::~ElementBndBoxTree() - { - for ( int i = 0; i < _elements.size(); ++i ) - if ( --_elements[i]->_refCount <= 0 ) - delete _elements[i]; - } - - //================================================================================ - /*! - * \brief Return the maximal box - */ - //================================================================================ - - Bnd_B3d* ElementBndBoxTree::buildRootBox() - { - Bnd_B3d* box = new Bnd_B3d; - for ( int i = 0; i < _elements.size(); ++i ) - box->Add( *_elements[i] ); - return box; - } - - //================================================================================ - /*! - * \brief Redistrubute element boxes among children - */ - //================================================================================ - - void ElementBndBoxTree::buildChildrenData() - { - for ( int i = 0; i < _elements.size(); ++i ) - { - for (int j = 0; j < 8; j++) - { - if ( !_elements[i]->IsOut( *myChildren[j]->getBox() )) - { - _elements[i]->_refCount++; - ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]); - } - } - _elements[i]->_refCount--; - } - _size = _elements.size(); - SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory - - for (int j = 0; j < 8; j++) - { - ElementBndBoxTree* child = static_cast( myChildren[j]); - if ( child->_elements.size() <= MaxNbElemsInLeaf ) - child->myIsLeaf = true; - - if ( child->_elements.capacity() - child->_elements.size() > 1000 ) - SMESHUtils::CompactVector( child->_elements ); - } - } - - //================================================================================ - /*! - * \brief Return elements which can include the point - */ - //================================================================================ - - void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point, - TIDSortedElemSet& foundElems) - { - if ( getBox()->IsOut( point.XYZ() )) - return; - - if ( isLeaf() ) - { - for ( int i = 0; i < _elements.size(); ++i ) - if ( !_elements[i]->IsOut( point.XYZ() )) - foundElems.insert( _elements[i]->_element ); - } - else - { - for (int i = 0; i < 8; i++) - ((ElementBndBoxTree*) myChildren[i])->getElementsNearPoint( point, foundElems ); - } - } - - //================================================================================ - /*! - * \brief Return elements which can be intersected by the line - */ - //================================================================================ - - void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line, - TIDSortedElemSet& foundElems) - { - if ( getBox()->IsOut( line )) - return; - - if ( isLeaf() ) - { - for ( int i = 0; i < _elements.size(); ++i ) - if ( !_elements[i]->IsOut( line )) - foundElems.insert( _elements[i]->_element ); - } - else - { - for (int i = 0; i < 8; i++) - ((ElementBndBoxTree*) myChildren[i])->getElementsNearLine( line, foundElems ); - } - } - - //================================================================================ - /*! - * \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 - */ - //================================================================================ - - ElementBndBoxTree::ElementBox::ElementBox(const SMDS_MeshElement* elem, double tolerance) - { - _element = elem; - _refCount = 1; - SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); - while ( nIt->more() ) - Add( SMESH_TNodeXYZ( nIt->next() )); - Enlarge( tolerance ); - } - -} // namespace - -//======================================================================= -/*! - * \brief Implementation of search for the elements by point and - * of classification of point in 2D mesh - */ -//======================================================================= - -struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher -{ - SMESHDS_Mesh* _mesh; - SMDS_ElemIteratorPtr _meshPartIt; - ElementBndBoxTree* _ebbTree; - SMESH_NodeSearcherImpl* _nodeSearcher; - SMDSAbs_ElementType _elementType; - double _tolerance; - bool _outerFacesFound; - set _outerFaces; // empty means "no internal faces at all" - - SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh, SMDS_ElemIteratorPtr elemIt=SMDS_ElemIteratorPtr()) - : _mesh(&mesh),_meshPartIt(elemIt),_ebbTree(0),_nodeSearcher(0),_tolerance(-1),_outerFacesFound(false) {} - ~SMESH_ElementSearcherImpl() - { - if ( _ebbTree ) delete _ebbTree; _ebbTree = 0; - if ( _nodeSearcher ) delete _nodeSearcher; _nodeSearcher = 0; - } - virtual int FindElementsByPoint(const gp_Pnt& point, - 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, - vector< const SMDS_MeshElement* >& foundElems); - double getTolerance(); - bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face, - const double tolerance, double & param); - void findOuterBoundary(const SMDS_MeshElement* anyOuterFace); - bool isOuterBoundary(const SMDS_MeshElement* face) const - { - return _outerFaces.empty() || _outerFaces.count(face); - } - struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState()) - { - const SMDS_MeshElement* _face; - gp_Vec _faceNorm; - bool _coincides; //!< the line lays in face plane - TInters(const SMDS_MeshElement* face, const gp_Vec& faceNorm, bool coinc=false) - : _face(face), _faceNorm( faceNorm ), _coincides( coinc ) {} - }; - struct TFaceLink //!< link and faces sharing it (used in findOuterBoundary()) - { - SMESH_TLink _link; - TIDSortedElemSet _faces; - TFaceLink( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshElement* face) - : _link( n1, n2 ), _faces( &face, &face + 1) {} - }; -}; - -ostream& operator<< (ostream& out, const SMESH_ElementSearcherImpl::TInters& i) -{ - return out << "TInters(face=" << ( i._face ? i._face->GetID() : 0) - << ", _coincides="<GetMeshInfo(); - - _tolerance = 0; - if ( _nodeSearcher && meshInfo.NbNodes() > 1 ) - { - double boxSize = _nodeSearcher->getTree()->maxSize(); - _tolerance = 1e-8 * boxSize/* / meshInfo.NbNodes()*/; - } - else if ( _ebbTree && meshInfo.NbElements() > 0 ) - { - double boxSize = _ebbTree->maxSize(); - _tolerance = 1e-8 * boxSize/* / meshInfo.NbElements()*/; - } - if ( _tolerance == 0 ) - { - // define tolerance by size of a most complex element - int complexType = SMDSAbs_Volume; - while ( complexType > SMDSAbs_All && - meshInfo.NbElements( SMDSAbs_ElementType( complexType )) < 1 ) - --complexType; - if ( complexType == SMDSAbs_All ) return 0; // empty mesh - double elemSize; - if ( complexType == int( SMDSAbs_Node )) - { - SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator(); - elemSize = 1; - if ( meshInfo.NbNodes() > 2 ) - elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() ); - } - else - { - SMDS_ElemIteratorPtr elemIt = - _mesh->elementsIterator( SMDSAbs_ElementType( complexType )); - const SMDS_MeshElement* elem = elemIt->next(); - SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); - SMESH_TNodeXYZ n1( cast2Node( nodeIt->next() )); - elemSize = 0; - while ( nodeIt->more() ) - { - double dist = n1.Distance( cast2Node( nodeIt->next() )); - elemSize = max( dist, elemSize ); - } - } - _tolerance = 1e-4 * elemSize; - } - } - return _tolerance; -} - -//================================================================================ -/*! - * \brief Find intersection of the line and an edge of face and return parameter on line - */ -//================================================================================ - -bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& line, - const SMDS_MeshElement* face, - const double tol, - double & param) -{ - int nbInts = 0; - param = 0; - - GeomAPI_ExtremaCurveCurve anExtCC; - Handle(Geom_Curve) lineCurve = new Geom_Line( line ); - - int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes(); - for ( int i = 0; i < nbNodes && nbInts < 2; ++i ) - { - GC_MakeSegment edge( SMESH_TNodeXYZ( face->GetNode( i )), - SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); - anExtCC.Init( lineCurve, edge); - if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol) - { - Quantity_Parameter pl, pe; - anExtCC.LowerDistanceParameters( pl, pe ); - param += pl; - if ( ++nbInts == 2 ) - break; - } - } - if ( nbInts > 0 ) param /= nbInts; - return nbInts > 0; -} -//================================================================================ -/*! - * \brief Find all faces belonging to the outer boundary of mesh - */ -//================================================================================ - -void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerFace) -{ - if ( _outerFacesFound ) return; - - // Collect all outer faces by passing from one outer face to another via their links - // and BTW find out if there are internal faces at all. - - // checked links and links where outer boundary meets internal one - set< SMESH_TLink > visitedLinks, seamLinks; - - // links to treat with already visited faces sharing them - list < TFaceLink > startLinks; - - // load startLinks with the first outerFace - startLinks.push_back( TFaceLink( outerFace->GetNode(0), outerFace->GetNode(1), outerFace)); - _outerFaces.insert( outerFace ); - - TIDSortedElemSet emptySet; - while ( !startLinks.empty() ) - { - const SMESH_TLink& link = startLinks.front()._link; - TIDSortedElemSet& faces = startLinks.front()._faces; - - outerFace = *faces.begin(); - // find other faces sharing the link - const SMDS_MeshElement* f; - while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces ))) - faces.insert( f ); - - // select another outer face among the found - const SMDS_MeshElement* outerFace2 = 0; - if ( faces.size() == 2 ) - { - outerFace2 = (outerFace == *faces.begin() ? *faces.rbegin() : *faces.begin()); - } - else if ( faces.size() > 2 ) - { - seamLinks.insert( link ); - - // link direction within the outerFace - gp_Vec n1n2( SMESH_TNodeXYZ( link.node1()), - SMESH_TNodeXYZ( link.node2())); - int i1 = outerFace->GetNodeIndex( link.node1() ); - int i2 = outerFace->GetNodeIndex( link.node2() ); - bool rev = ( abs(i2-i1) == 1 ? i1 > i2 : i2 > i1 ); - if ( rev ) n1n2.Reverse(); - // outerFace normal - gp_XYZ ofNorm, fNorm; - if ( SMESH_Algo::FaceNormal( outerFace, ofNorm, /*normalized=*/false )) - { - // direction from the link inside outerFace - gp_Vec dirInOF = gp_Vec( ofNorm ) ^ n1n2; - // sort all other faces by angle with the dirInOF - map< double, const SMDS_MeshElement* > angle2Face; - set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin(); - for ( ; face != faces.end(); ++face ) - { - if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false )) - continue; - gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2; - double angle = dirInOF.AngleWithRef( dirInF, n1n2 ); - if ( angle < 0 ) angle += 2. * M_PI; - angle2Face.insert( make_pair( angle, *face )); - } - if ( !angle2Face.empty() ) - outerFace2 = angle2Face.begin()->second; - } - } - // store the found outer face and add its links to continue seaching from - if ( outerFace2 ) - { - _outerFaces.insert( outerFace ); - int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 ); - for ( int i = 0; i < nbNodes; ++i ) - { - SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes)); - if ( visitedLinks.insert( link2 ).second ) - startLinks.push_back( TFaceLink( link2.node1(), link2.node2(), outerFace2 )); - } - } - startLinks.pop_front(); - } - _outerFacesFound = true; - - if ( !seamLinks.empty() ) - { - // There are internal boundaries touching the outher one, - // find all faces of internal boundaries in order to find - // faces of boundaries of holes, if any. - - } - else - { - _outerFaces.clear(); - } -} - -//======================================================================= -/*! - * \brief Find elements of given type where the given point is IN or ON. - * Returns nb of found elements and elements them-selves. - * - * 'ALL' type means elements of any type excluding nodes, balls and 0D elements - */ -//======================================================================= - -int SMESH_ElementSearcherImpl:: -FindElementsByPoint(const gp_Pnt& point, - SMDSAbs_ElementType type, - vector< const SMDS_MeshElement* >& foundElements) -{ - foundElements.clear(); - - double tolerance = getTolerance(); - - // ================================================================================= - if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball) - { - if ( !_nodeSearcher ) - _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh ); - - const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point ); - if ( !closeNode ) return foundElements.size(); - - if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance ) - return foundElements.size(); // to far from any node - - if ( type == SMDSAbs_Node ) - { - foundElements.push_back( closeNode ); - } - else - { - SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type ); - while ( elemIt->more() ) - foundElements.push_back( elemIt->next() ); - } - } - // ================================================================================= - else // elements more complex than 0D - { - if ( !_ebbTree || _elementType != type ) - { - if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance ); - } - TIDSortedElemSet suspectElems; - _ebbTree->getElementsNearPoint( point, suspectElems ); - TIDSortedElemSet::iterator elem = suspectElems.begin(); - for ( ; elem != suspectElems.end(); ++elem ) - 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::max(); - multimap< double, const SMDS_MeshElement* > dist2face; - TIDSortedElemSet::iterator elem = suspectElems.begin(); - for ( ; elem != suspectElems.end(); ++elem ) - { - double dist = SMESH_MeshEditor::GetDistance( dynamic_cast(*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( " <GetID() << " DIST " << minDist << endl; - } - } - else - { - // NOT IMPLEMENTED SO FAR - } - return closestElem; -} - - -//================================================================================ -/*! - * \brief Classify the given point in the closed 2D mesh - */ -//================================================================================ - -TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) -{ - double tolerance = getTolerance(); - if ( !_ebbTree || _elementType != SMDSAbs_Face ) - { - if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face, _meshPartIt ); - } - // Algo: analyse transition of a line starting at the point through mesh boundary; - // try three lines parallel to axis of the coordinate system and perform rough - // analysis. If solution is not clear perform thorough analysis. - - const int nbAxes = 3; - gp_Dir axisDir[ nbAxes ] = { gp::DX(), gp::DY(), gp::DZ() }; - map< double, TInters > paramOnLine2TInters[ nbAxes ]; - list< TInters > tangentInters[ nbAxes ]; // of faces whose plane includes the line - multimap< int, int > nbInt2Axis; // to find the simplest case - for ( int axis = 0; axis < nbAxes; ++axis ) - { - gp_Ax1 lineAxis( point, axisDir[axis]); - gp_Lin line ( lineAxis ); - - TIDSortedElemSet suspectFaces; // faces possibly intersecting the line - _ebbTree->getElementsNearLine( lineAxis, suspectFaces ); - - // Intersect faces with the line - - map< double, TInters > & u2inters = paramOnLine2TInters[ axis ]; - TIDSortedElemSet::iterator face = suspectFaces.begin(); - for ( ; face != suspectFaces.end(); ++face ) - { - // get face plane - gp_XYZ fNorm; - if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue; - gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm ); - - // perform intersection - IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane )); - if ( !intersection.IsDone() ) - continue; - if ( intersection.IsInQuadric() ) - { - tangentInters[ axis ].push_back( TInters( *face, fNorm, true )); - } - else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 ) - { - gp_Pnt intersectionPoint = intersection.Point(1); - if ( !SMESH_MeshEditor::IsOut( *face, intersectionPoint, tolerance )) - u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm ))); - } - } - // Analyse intersections roughly - - int nbInter = u2inters.size(); - if ( nbInter == 0 ) - return TopAbs_OUT; - - double f = u2inters.begin()->first, l = u2inters.rbegin()->first; - if ( nbInter == 1 ) // not closed mesh - return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN; - - if ( fabs( f ) < tolerance || fabs( l ) < tolerance ) - return TopAbs_ON; - - if ( (f<0) == (l<0) ) - return TopAbs_OUT; - - int nbIntBeforePoint = std::distance( u2inters.begin(), u2inters.lower_bound(0)); - int nbIntAfterPoint = nbInter - nbIntBeforePoint; - if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 ) - return TopAbs_IN; - - nbInt2Axis.insert( make_pair( min( nbIntBeforePoint, nbIntAfterPoint ), axis )); - - if ( _outerFacesFound ) break; // pass to thorough analysis - - } // three attempts - loop on CS axes - - // Analyse intersections thoroughly. - // We make two loops maximum, on the first one we only exclude touching intersections, - // on the second, if situation is still unclear, we gather and use information on - // position of faces (internal or outer). If faces position is already gathered, - // we make the second loop right away. - - for ( int hasPositionInfo = _outerFacesFound; hasPositionInfo < 2; ++hasPositionInfo ) - { - multimap< int, int >::const_iterator nb_axis = nbInt2Axis.begin(); - for ( ; nb_axis != nbInt2Axis.end(); ++nb_axis ) - { - int axis = nb_axis->second; - map< double, TInters > & u2inters = paramOnLine2TInters[ axis ]; - - gp_Ax1 lineAxis( point, axisDir[axis]); - gp_Lin line ( lineAxis ); - - // add tangent intersections to u2inters - double param; - list< TInters >::const_iterator tgtInt = tangentInters[ axis ].begin(); - for ( ; tgtInt != tangentInters[ axis ].end(); ++tgtInt ) - if ( getIntersParamOnLine( line, tgtInt->_face, tolerance, param )) - u2inters.insert(make_pair( param, *tgtInt )); - tangentInters[ axis ].clear(); - - // Count intersections before and after the point excluding touching ones. - // If hasPositionInfo we count intersections of outer boundary only - - int nbIntBeforePoint = 0, nbIntAfterPoint = 0; - double f = numeric_limits::max(), l = -numeric_limits::max(); - map< double, TInters >::iterator u_int1 = u2inters.begin(), u_int2 = u_int1; - bool ok = ! u_int1->second._coincides; - while ( ok && u_int1 != u2inters.end() ) - { - double u = u_int1->first; - bool touchingInt = false; - if ( ++u_int2 != u2inters.end() ) - { - // skip intersections at the same point (if the line passes through edge or node) - int nbSamePnt = 0; - while ( u_int2 != u2inters.end() && fabs( u_int2->first - u ) < tolerance ) - { - ++nbSamePnt; - ++u_int2; - } - - // skip tangent intersections - int nbTgt = 0; - const SMDS_MeshElement* prevFace = u_int1->second._face; - while ( ok && u_int2->second._coincides ) - { - if ( SMESH_Algo::GetCommonNodes(prevFace , u_int2->second._face).empty() ) - ok = false; - else - { - nbTgt++; - u_int2++; - ok = ( u_int2 != u2inters.end() ); - } - } - if ( !ok ) break; - - // skip intersections at the same point after tangent intersections - if ( nbTgt > 0 ) - { - double u2 = u_int2->first; - ++u_int2; - while ( u_int2 != u2inters.end() && fabs( u_int2->first - u2 ) < tolerance ) - { - ++nbSamePnt; - ++u_int2; - } - } - // decide if we skipped a touching intersection - if ( nbSamePnt + nbTgt > 0 ) - { - double minDot = numeric_limits::max(), maxDot = -numeric_limits::max(); - map< double, TInters >::iterator u_int = u_int1; - for ( ; u_int != u_int2; ++u_int ) - { - if ( u_int->second._coincides ) continue; - double dot = u_int->second._faceNorm * line.Direction(); - if ( dot > maxDot ) maxDot = dot; - if ( dot < minDot ) minDot = dot; - } - touchingInt = ( minDot*maxDot < 0 ); - } - } - if ( !touchingInt ) - { - if ( !hasPositionInfo || isOuterBoundary( u_int1->second._face )) - { - if ( u < 0 ) - ++nbIntBeforePoint; - else - ++nbIntAfterPoint; - } - if ( u < f ) f = u; - if ( u > l ) l = u; - } - - u_int1 = u_int2; // to next intersection - - } // loop on intersections with one line - - if ( ok ) - { - if ( fabs( f ) < tolerance || fabs( l ) < tolerance ) - return TopAbs_ON; - - if ( nbIntBeforePoint == 0 || nbIntAfterPoint == 0) - return TopAbs_OUT; - - if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh - return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN; - - if ( nbIntBeforePoint == 1 || nbIntAfterPoint == 1 ) - return TopAbs_IN; - - if ( (f<0) == (l<0) ) - return TopAbs_OUT; - - if ( hasPositionInfo ) - return nbIntBeforePoint % 2 ? TopAbs_IN : TopAbs_OUT; - } - } // loop on intersections of the tree lines - thorough analysis - - if ( !hasPositionInfo ) - { - // gather info on faces position - is face in the outer boundary or not - map< double, TInters > & u2inters = paramOnLine2TInters[ 0 ]; - findOuterBoundary( u2inters.begin()->second._face ); - } - - } // two attempts - with and w/o faces position info in the mesh - - return TopAbs_UNKNOWN; -} - -//======================================================================= -/*! - * \brief Return elements possibly intersecting the line - */ -//======================================================================= - -void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& line, - SMDSAbs_ElementType type, - vector< const SMDS_MeshElement* >& foundElems) -{ - if ( !_ebbTree || _elementType != type ) - { - if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt ); - } - TIDSortedElemSet suspectFaces; // elements possibly intersecting the line - _ebbTree->getElementsNearLine( line, suspectFaces ); - foundElems.assign( suspectFaces.begin(), suspectFaces.end()); -} - -//======================================================================= -/*! - * \brief Return SMESH_ElementSearcher - */ -//======================================================================= - -SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher() -{ - return new SMESH_ElementSearcherImpl( *GetMeshDS() ); -} - -//======================================================================= -/*! - * \brief Return SMESH_ElementSearcher acting on a sub-set of elements - */ -//======================================================================= - -SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher(SMDS_ElemIteratorPtr elemIt) -{ - return new SMESH_ElementSearcherImpl( *GetMeshDS(), elemIt ); -} - -//======================================================================= -/*! - * \brief Return true if the point is IN or ON of the element - */ -//======================================================================= - -bool SMESH_MeshEditor::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol ) -{ - if ( element->GetType() == SMDSAbs_Volume) - { - return SMDS_VolumeTool( element ).IsOut( point.X(), point.Y(), point.Z(), tol ); - } - - // get ordered nodes - - vector< gp_XYZ > xyz; - vector nodeList; - - SMDS_ElemIteratorPtr nodeIt = element->nodesIterator(); - if ( element->IsQuadratic() ) { - if (const SMDS_VtkFace* f=dynamic_cast(element)) - nodeIt = f->interlacedNodesElemIterator(); - else if (const SMDS_VtkEdge* e =dynamic_cast(element)) - nodeIt = e->interlacedNodesElemIterator(); - } - while ( nodeIt->more() ) - { - const SMDS_MeshNode* node = cast2Node( nodeIt->next() ); - xyz.push_back( SMESH_TNodeXYZ(node) ); - nodeList.push_back(node); - } - - int i, nbNodes = element->NbNodes(); - - if ( element->GetType() == SMDSAbs_Face ) // -------------------------------------------------- - { - // compute face normal - gp_Vec faceNorm(0,0,0); - xyz.push_back( xyz.front() ); - nodeList.push_back( nodeList.front() ); - for ( i = 0; i < nbNodes; ++i ) - { - gp_Vec edge1( xyz[i+1], xyz[i]); - gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] ); - faceNorm += edge1 ^ edge2; - } - double normSize = faceNorm.Magnitude(); - if ( normSize <= tol ) - { - // degenerated face: point is out if it is out of all face edges - for ( i = 0; i < nbNodes; ++i ) + // there must be a top element + const SMDS_MeshElement* topElem = 0; + if ( isNodes ) { - SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] ); - if ( !IsOut( &edge, point, tol )) - return false; + topElem = resultElems.back(); + resultElems.pop_back(); } - return true; - } - faceNorm /= normSize; - - // check if the point lays on face plane - gp_Vec n2p( xyz[0], point ); - if ( fabs( n2p * faceNorm ) > tol ) - return true; // not on face plane - - // check if point is out of face boundary: - // define it by closest transition of a ray point->infinity through face boundary - // on the face plane. - // First, find normal of a plane perpendicular to face plane, to be used as a cutting tool - // to find intersections of the ray with the boundary. - gp_Vec ray = n2p; - gp_Vec plnNorm = ray ^ faceNorm; - normSize = plnNorm.Magnitude(); - if ( normSize <= tol ) return false; // point coincides with the first node - plnNorm /= normSize; - // for each node of the face, compute its signed distance to the plane - vector dist( nbNodes + 1); - for ( i = 0; i < nbNodes; ++i ) - { - gp_Vec n2p( xyz[i], point ); - dist[i] = n2p * plnNorm; - } - dist.back() = dist.front(); - // find the closest intersection - int iClosest = -1; - double rClosest, distClosest = 1e100;; - gp_Pnt pClosest; - for ( i = 0; i < nbNodes; ++i ) - { - double r; - if ( fabs( dist[i]) < tol ) - r = 0.; - else if ( fabs( dist[i+1]) < tol ) - r = 1.; - else if ( dist[i] * dist[i+1] < 0 ) - r = dist[i] / ( dist[i] - dist[i+1] ); else - continue; // no intersection - gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r; - gp_Vec p2int ( point, pInt); - if ( p2int * ray > -tol ) // right half-space { - double intDist = p2int.SquareMagnitude(); - if ( intDist < distClosest ) + list< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin(); + for ( ; resElemIt != resultElems.rend() ; ++resElemIt ) + if ( (*resElemIt)->GetType() == sourceElem->GetType() ) + { + topElem = *resElemIt; + resultElems.erase( --(resElemIt.base()) ); // erase *resElemIt + break; + } + } + + // add resultElems to groups originted from ones the sourceElem belongs to + list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end(); + for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew ) + { + SMESHDS_GroupBase* oldGroup = gOldNew->get<0>(); + if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup { - iClosest = i; - rClosest = r; - pClosest = pInt; - distClosest = intDist; + // fill in a new group + SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup(); + list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt; + for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt ) + newGroup.Add( *resElemIt ); + + // fill a "top" group + if ( topElem ) + { + SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup(); + newTopGroup.Add( topElem ); + } } } - } - if ( iClosest < 0 ) - return true; // no intesections - out - - // analyse transition - gp_Vec edge( xyz[iClosest], xyz[iClosest+1] ); - gp_Vec edgeNorm = -( edge ^ faceNorm ); // normal to intersected edge pointing out of face - gp_Vec p2int ( point, pClosest ); - bool out = (edgeNorm * p2int) < -tol; - if ( rClosest > 0. && rClosest < 1. ) // not node intersection - return out; + } // loop on created elements + }// loop on nodes and elements - // ray pass through a face node; analyze transition through an adjacent edge - gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ]; - gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ]; - gp_Vec edgeAdjacent( p1, p2 ); - gp_Vec edgeNorm2 = -( edgeAdjacent ^ faceNorm ); - bool out2 = (edgeNorm2 * p2int) < -tol; + // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups - bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0; - return covexCorner ? (out || out2) : (out && out2); - } - if ( element->GetType() == SMDSAbs_Edge ) // -------------------------------------------------- + list topGrouIds; + for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i ) { - // point is out of edge if it is NOT ON any straight part of edge - // (we consider quadratic edge as being composed of two straight parts) - for ( i = 1; i < nbNodes; ++i ) + SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->get<0>(); + SMESHDS_Group* newGroups[2] = { orderedOldNewGroups[i]->get<1>(), + orderedOldNewGroups[i]->get<2>() }; + const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty(); + for ( int is2nd = 0; is2nd < 2; ++is2nd ) { - gp_Vec edge( xyz[i-1], xyz[i]); - gp_Vec n1p ( xyz[i-1], point); - double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude(); - if ( dist > tol ) - continue; - gp_Vec n2p( xyz[i], point ); - if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol ) - continue; - return false; // point is ON this part - } - return true; - } - // Node or 0D element ------------------------------------------------------------------------- - { - gp_Vec n2p ( xyz[0], point ); - return n2p.Magnitude() <= tol; - } - return true; -} + SMESHDS_Group* newGroupDS = newGroups[ is2nd ]; + if ( newGroupDS->IsEmpty() ) + { + mesh->GetMeshDS()->RemoveGroup( newGroupDS ); + } + else + { + // set group type + newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() ); -//======================================================================= + // make a name + const bool isTop = ( nbNewGroups == 2 && + newGroupDS->GetType() == oldGroupDS->GetType() && + is2nd ); -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 + string name = oldGroupDS->GetStoreName(); + if ( !targetMesh ) { + string suffix = ( isTop ? "top": postfix.c_str() ); + name += "_"; + name += suffix; + int nb = 1; + while ( !groupNames.insert( name ).second ) // name exists + name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << suffix << "_" << nb++; + } + else if ( isTop ) { + name += "_top"; + } + newGroupDS->SetStoreName( name.c_str() ); - 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; + // make a SMESH_Groups + mesh->AddGroup( newGroupDS ); + if ( isTop ) + topGrouIds.push_back( newGroupDS->GetID() ); + else + newGroupIDs->push_back( newGroupDS->GetID() ); + } } - }; - - //================================================================================ - /*! - * \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 ); } + newGroupIDs->splice( newGroupIDs->end(), topGrouIds ); + + return newGroupIDs; } -//======================================================================= +//================================================================================ /*! - * \brief Return minimal distance from a point to a face - * - * Currently we ignore non-planarity and 2nd order of face + * \brief Return list of group of nodes close to each other within theTolerance + * Search among theNodes or in the whole mesh if theNodes is empty using + * an Octree algorithm */ -//======================================================================= +//================================================================================ -double SMESH_MeshEditor::GetDistance( const SMDS_MeshFace* face, - const gp_Pnt& point ) +void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes, + const double theTolerance, + TListOfListOfNodes & theGroupsOfNodes) { - double badDistance = -1; - if ( !face ) return badDistance; - - // coordinates of nodes (medium nodes, if any, ignored) - typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; - vector 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::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 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() ); + myLastCreatedElems.Clear(); + myLastCreatedNodes.Clear(); - // 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 ); + if ( theNodes.empty() ) + { // get all nodes in the mesh + SMDS_NodeIteratorPtr nIt = GetMeshDS()->nodesIterator(/*idInceasingOrder=*/true); + while ( nIt->more() ) + theNodes.insert( theNodes.end(),nIt->next()); } - // 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; + SMESH_OctreeNode::FindCoincidentNodes ( theNodes, &theGroupsOfNodes, theTolerance); } //======================================================================= //function : SimplifyFace //purpose : //======================================================================= -int SMESH_MeshEditor::SimplifyFace (const vector faceNodes, - vector& poly_nodes, - vector& quantities) const + +int SMESH_MeshEditor::SimplifyFace (const vector& faceNodes, + vector& poly_nodes, + vector& quantities) const { int nbNodes = faceNodes.size(); @@ -8414,75 +7131,6 @@ void SMESH_MeshEditor::MergeEqualElements() MergeElements(aGroupsOfElementsID); } -//======================================================================= -//function : FindFaceInSet -//purpose : Return a face having linked nodes n1 and n2 and which is -// - not in avoidSet, -// - in elemSet provided that !elemSet.empty() -// i1 and i2 optionally returns indices of n1 and n2 -//======================================================================= - -const SMDS_MeshElement* -SMESH_MeshEditor::FindFaceInSet(const SMDS_MeshNode* n1, - const SMDS_MeshNode* n2, - const TIDSortedElemSet& elemSet, - const TIDSortedElemSet& avoidSet, - int* n1ind, - int* n2ind) - -{ - int i1, i2; - const SMDS_MeshElement* face = 0; - - SMDS_ElemIteratorPtr invElemIt = n1->GetInverseElementIterator(SMDSAbs_Face); - //MESSAGE("n1->GetInverseElementIterator(SMDSAbs_Face) " << invElemIt); - while ( invElemIt->more() && !face ) // loop on inverse faces of n1 - { - //MESSAGE("in while ( invElemIt->more() && !face )"); - const SMDS_MeshElement* elem = invElemIt->next(); - if (avoidSet.count( elem )) - continue; - if ( !elemSet.empty() && !elemSet.count( elem )) - continue; - // index of n1 - i1 = elem->GetNodeIndex( n1 ); - // find a n2 linked to n1 - int nbN = elem->IsQuadratic() ? elem->NbNodes()/2 : elem->NbNodes(); - for ( int di = -1; di < 2 && !face; di += 2 ) - { - i2 = (i1+di+nbN) % nbN; - if ( elem->GetNode( i2 ) == n2 ) - face = elem; - } - if ( !face && elem->IsQuadratic()) - { - // analysis for quadratic elements using all nodes - const SMDS_VtkFace* F = - dynamic_cast(elem); - if (!F) throw SALOME_Exception(LOCALIZED("not an SMDS_VtkFace")); - // use special nodes iterator - SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator(); - const SMDS_MeshNode* prevN = cast2Node( anIter->next() ); - for ( i1 = -1, i2 = 0; anIter->more() && !face; i1++, i2++ ) - { - const SMDS_MeshNode* n = cast2Node( anIter->next() ); - if ( n1 == prevN && n2 == n ) - { - face = elem; - } - else if ( n2 == prevN && n1 == n ) - { - face = elem; swap( i1, i2 ); - } - prevN = n; - } - } - } - if ( n1ind ) *n1ind = i1; - if ( n2ind ) *n2ind = i2; - return face; -} - //======================================================================= //function : findAdjacentFace //purpose : @@ -8495,7 +7143,7 @@ static const SMDS_MeshElement* findAdjacentFace(const SMDS_MeshNode* n1, TIDSortedElemSet elemSet, avoidSet; if ( elem ) avoidSet.insert ( elem ); - return SMESH_MeshEditor::FindFaceInSet( n1, n2, elemSet, avoidSet ); + return SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet ); } //======================================================================= @@ -8611,7 +7259,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode* theFirst cNL = & contNodes[ contNodes[0].empty() ? 0 : 1 ]; cFL = & contFaces[ contFaces[0].empty() ? 0 : 1 ]; // find one more free border - if ( ! FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) { + if ( ! SMESH_MeshEditor::FindFreeBorder( nStart, *nStartIt, theLastNode, *cNL, *cFL )) { cNL->clear(); cFL->clear(); } @@ -9595,24 +8243,47 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, const SMDS_MeshElement* elem = ElemItr->next(); if( !elem ) continue; + // analyse a necessity of conversion + const SMDSAbs_ElementType aType = elem->GetType(); + if ( aType < SMDSAbs_Edge || aType > SMDSAbs_Volume ) + continue; const SMDSAbs_EntityType aGeomType = elem->GetEntityType(); + bool hasCentralNodes = false; if ( elem->IsQuadratic() ) { bool alreadyOK; switch ( aGeomType ) { + case SMDSEntity_Quad_Triangle: case SMDSEntity_Quad_Quadrangle: - case SMDSEntity_Quad_Hexa: alreadyOK = !theHelper.GetIsBiQuadratic(); break; + case SMDSEntity_Quad_Hexa: + alreadyOK = !theHelper.GetIsBiQuadratic(); break; + + case SMDSEntity_BiQuad_Triangle: case SMDSEntity_BiQuad_Quadrangle: - case SMDSEntity_TriQuad_Hexa: alreadyOK = theHelper.GetIsBiQuadratic(); break; - default: alreadyOK = true; - } - if ( alreadyOK ) continue; + case SMDSEntity_TriQuad_Hexa: + alreadyOK = theHelper.GetIsBiQuadratic(); + hasCentralNodes = true; + break; + default: + alreadyOK = true; + } + // take into account already present modium nodes + switch ( aType ) { + case SMDSAbs_Volume: + theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break; + case SMDSAbs_Face: + theHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( elem )); break; + case SMDSAbs_Edge: + theHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( elem )); break; + default:; + } + if ( alreadyOK ) + continue; } // get elem data needed to re-create it // - const int id = elem->GetID(); - const int nbNodes = elem->NbCornerNodes(); - const SMDSAbs_ElementType aType = elem->GetType(); + const int id = elem->GetID(); + const int nbNodes = elem->NbCornerNodes(); nodes.assign(elem->begin_nodes(), elem->end_nodes()); if ( aGeomType == SMDSEntity_Polyhedra ) nbNodeInFaces = static_cast( elem )->GetQuantities(); @@ -9622,6 +8293,12 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, // remove a linear element GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false); + // remove central nodes of biquadratic elements (biquad->quad convertion) + if ( hasCentralNodes ) + for ( size_t i = nbNodes * 2; i < nodes.size(); ++i ) + if ( nodes[i]->NbInverseElements() == 0 ) + GetMeshDS()->RemoveFreeNode( nodes[i], theSm, /*fromGroups=*/true ); + const SMDS_MeshElement* NewElem = 0; switch( aType ) @@ -9643,7 +8320,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, break; default: NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d); - continue; } break; } @@ -9676,7 +8352,7 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, continue; } ReplaceElemInGroups( elem, NewElem, GetMeshDS()); - if( NewElem ) + if( NewElem && NewElem->getshapeId() < 1 ) theSm->AddElement( NewElem ); } return nbElem; @@ -9696,6 +8372,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT aHelper.SetIsBiQuadratic( theToBiQuad ); aHelper.SetElementsOnShape(true); + // convert elements assigned to sub-meshes int nbCheckedElems = 0; if ( myMesh->HasShapeToMesh() ) { @@ -9711,19 +8388,22 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT } } } + + // convert elements NOT assigned to sub-meshes int totalNbElems = meshDS->NbEdges() + meshDS->NbFaces() + meshDS->NbVolumes(); - if ( nbCheckedElems < totalNbElems ) // not all elements are in submeshes + if ( nbCheckedElems < totalNbElems ) // not all elements are in sub-meshes { aHelper.SetElementsOnShape(false); SMESHDS_SubMesh *smDS = 0; + + // convert edges SMDS_EdgeIteratorPtr aEdgeItr = meshDS->edgesIterator(); - while(aEdgeItr->more()) + while( aEdgeItr->more() ) { const SMDS_MeshEdge* edge = aEdgeItr->next(); - if(edge && !edge->IsQuadratic()) + if ( !edge->IsQuadratic() ) { - int id = edge->GetID(); - //MESSAGE("edge->GetID() " << id); + int id = edge->GetID(); const SMDS_MeshNode* n1 = edge->GetNode(0); const SMDS_MeshNode* n2 = edge->GetNode(1); @@ -9732,16 +8412,36 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT const SMDS_MeshEdge* NewEdge = aHelper.AddEdge(n1, n2, id, theForce3d); ReplaceElemInGroups( edge, NewEdge, GetMeshDS()); } + else + { + aHelper.AddTLinks( static_cast< const SMDS_MeshEdge* >( edge )); + } } + + // convert faces SMDS_FaceIteratorPtr aFaceItr = meshDS->facesIterator(); - while(aFaceItr->more()) + while( aFaceItr->more() ) { const SMDS_MeshFace* face = aFaceItr->next(); if ( !face ) continue; const SMDSAbs_EntityType type = face->GetEntityType(); - if (( theToBiQuad && type == SMDSEntity_BiQuad_Quadrangle ) || - ( !theToBiQuad && type == SMDSEntity_Quad_Quadrangle )) + bool alreadyOK; + switch( type ) + { + case SMDSEntity_Quad_Triangle: + case SMDSEntity_Quad_Quadrangle: + alreadyOK = !theToBiQuad; + aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face )); + break; + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_BiQuad_Quadrangle: + alreadyOK = theToBiQuad; + aHelper.AddTLinks( static_cast< const SMDS_MeshFace* >( face )); + break; + default: alreadyOK = false; + } + if ( alreadyOK ) continue; const int id = face->GetID(); @@ -9753,28 +8453,42 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT switch( type ) { case SMDSEntity_Triangle: + case SMDSEntity_Quad_Triangle: + case SMDSEntity_BiQuad_Triangle: NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], id, theForce3d); + if ( nodes.size() == 7 && nodes[6]->NbInverseElements() == 0 ) // rm a central node + GetMeshDS()->RemoveFreeNode( nodes[6], /*sm=*/0, /*fromGroups=*/true ); break; + case SMDSEntity_Quadrangle: + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_BiQuad_Quadrangle: NewFace = aHelper.AddFace(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d); + if ( nodes.size() == 9 && nodes[8]->NbInverseElements() == 0 ) // rm a central node + GetMeshDS()->RemoveFreeNode( nodes[8], /*sm=*/0, /*fromGroups=*/true ); break; - default: + + default:; NewFace = aHelper.AddPolygonalFace(nodes, id, theForce3d); } ReplaceElemInGroups( face, NewFace, GetMeshDS()); } + + // convert volumes vector nbNodeInFaces; SMDS_VolumeIteratorPtr aVolumeItr = meshDS->volumesIterator(); while(aVolumeItr->more()) { const SMDS_MeshVolume* volume = aVolumeItr->next(); - if(!volume || volume->IsQuadratic() ) continue; + if ( !volume ) continue; const SMDSAbs_EntityType type = volume->GetEntityType(); if (( theToBiQuad && type == SMDSEntity_TriQuad_Hexa ) || ( !theToBiQuad && type == SMDSEntity_Quad_Hexa )) + { + aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume )); continue; - + } const int id = volume->GetID(); vector nodes (volume->begin_nodes(), volume->end_nodes()); if ( type == SMDSEntity_Polyhedra ) @@ -9795,6 +8509,9 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT case SMDSEntity_TriQuad_Hexa: NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d); + for ( size_t i = 20; i < nodes.size(); ++i ) // rm central nodes + if ( nodes[i]->NbInverseElements() == 0 ) + GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true ); break; case SMDSEntity_Pyramid: NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2], @@ -9814,8 +8531,9 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT if ( !theForce3d ) { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion - aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh - aHelper.FixQuadraticElements(myError); + // aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh + // aHelper.FixQuadraticElements(myError); + SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError); } } @@ -9853,36 +8571,36 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator(); while ( invIt->more() ) { - const SMDS_MeshElement* e = invIt->next(); + const SMDS_MeshElement* e = invIt->next(); + const SMDSAbs_ElementType type = e->GetType(); if ( e->IsQuadratic() ) { + quadAdjacentElems[ type ].insert( e ); + bool alreadyOK; switch ( e->GetEntityType() ) { + case SMDSEntity_Quad_Triangle: case SMDSEntity_Quad_Quadrangle: case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break; + case SMDSEntity_BiQuad_Triangle: case SMDSEntity_BiQuad_Quadrangle: case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break; default: alreadyOK = true; } if ( alreadyOK ) - { - quadAdjacentElems[ e->GetType() ].insert( e ); continue; - } - } - if ( e->GetType() >= elemType ) - { - continue; // same type of more complex linear element } + if ( type >= elemType ) + continue; // same type or more complex linear element - if ( !checkedAdjacentElems[ e->GetType() ].insert( e ).second ) + if ( !checkedAdjacentElems[ type ].insert( e ).second ) continue; // e is already checked // check nodes bool allIn = true; - SMDS_ElemIteratorPtr nodeIt = e->nodesIterator(); + SMDS_NodeIteratorPtr nodeIt = e->nodeIterator(); while ( nodeIt->more() && allIn ) - allIn = allNodes.count( cast2Node( nodeIt->next() )); + allIn = allNodes.count( nodeIt->next() ); if ( allIn ) theElements.insert(e ); } @@ -9920,27 +8638,38 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt ) { const SMDS_MeshElement* elem = *eIt; - if( elem->NbNodes() < 2 || elem->IsPoly() ) - continue; - if ( elem->IsQuadratic() ) - { - bool alreadyOK; - switch ( elem->GetEntityType() ) { - case SMDSEntity_Quad_Quadrangle: - case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break; - case SMDSEntity_BiQuad_Quadrangle: - case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break; - default: alreadyOK = true; - } - if ( alreadyOK ) continue; - } + bool alreadyOK; + int nbCentralNodes = 0; + switch ( elem->GetEntityType() ) { + // linear convertible + case SMDSEntity_Edge: + case SMDSEntity_Triangle: + case SMDSEntity_Quadrangle: + case SMDSEntity_Tetra: + case SMDSEntity_Pyramid: + case SMDSEntity_Hexa: + case SMDSEntity_Penta: alreadyOK = false; nbCentralNodes = 0; break; + // quadratic that can become bi-quadratic + case SMDSEntity_Quad_Triangle: + case SMDSEntity_Quad_Quadrangle: + case SMDSEntity_Quad_Hexa: alreadyOK =!theToBiQuad; nbCentralNodes = 0; break; + // bi-quadratic + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_BiQuad_Quadrangle: alreadyOK = theToBiQuad; nbCentralNodes = 1; break; + case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; nbCentralNodes = 7; break; + // the rest + default: alreadyOK = true; + } + if ( alreadyOK ) continue; const SMDSAbs_ElementType type = elem->GetType(); const int id = elem->GetID(); const int nbNodes = elem->NbCornerNodes(); vector nodes ( elem->begin_nodes(), elem->end_nodes()); + helper.SetSubShape( elem->getshapeId() ); + if ( !smDS || !smDS->Contains( elem )) smDS = meshDS->MeshElements( elem->getshapeId() ); meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false); @@ -9977,12 +8706,19 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, ReplaceElemInGroups( elem, newElem, meshDS); if( newElem && smDS ) smDS->AddElement( newElem ); - } - if ( !theForce3d && !getenv("NO_FixQuadraticElements")) + // remove central nodes + for ( size_t i = nodes.size() - nbCentralNodes; i < nodes.size(); ++i ) + if ( nodes[i]->NbInverseElements() == 0 ) + meshDS->RemoveFreeNode( nodes[i], smDS, /*fromGroups=*/true ); + + } // loop on theElements + + if ( !theForce3d ) { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion - helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh - helper.FixQuadraticElements( myError ); + // helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh + // helper.FixQuadraticElements( myError ); + SMESH_MesherHelper( *myMesh ).FixQuadraticElements(myError); } } @@ -10103,8 +8839,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements) } // replace given elements by linear ones - typedef SMDS_SetIterator TSetIterator; - SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() )); + SMDS_ElemIteratorPtr elemIt = elemSetIterator( theElements ); removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 ); // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs @@ -10141,7 +8876,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements) const SMDS_MeshElement* eComplex = invIt2->next(); if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes)) { - int nbCommonNodes = SMESH_Algo::GetCommonNodes( e, eComplex ).size(); + int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e, eComplex ).size(); if ( nbCommonNodes == e->NbNodes()) { complexFound = true; @@ -10156,8 +8891,7 @@ void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements) } } } - elemIt = SMDS_ElemIteratorPtr - (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() )); + elemIt = elemSetIterator( moreElemsToConvert ); removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 ); } @@ -10506,8 +9240,10 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1, //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] ); + face[ iSide ] = SMESH_MeshAlgos::FindFaceInSet( n1, n2, + *faceSetPtr[ iSide ], avoidSet, + &iLinkNode[iSide][0], + &iLinkNode[iSide][1] ); if ( face[ iSide ]) { //cout << " F " << face[ iSide]->GetID() <& theSide1, return SEW_OK; } +//================================================================================ +/*! + * \brief Create elements equal (on same nodes) to given ones + * \param [in] theElements - a set of elems to duplicate. If it is empty, all + * elements of the uppest dimension are duplicated. + */ +//================================================================================ + +void SMESH_MeshEditor::DoubleElements( const TIDSortedElemSet& theElements ) +{ + CrearLastCreated(); + SMESHDS_Mesh* mesh = GetMeshDS(); + + // get an element type and an iterator over elements + + SMDSAbs_ElementType type; + SMDS_ElemIteratorPtr elemIt; + vector< const SMDS_MeshElement* > allElems; + if ( theElements.empty() ) + { + if ( mesh->NbNodes() == 0 ) + return; + // get most complex type + SMDSAbs_ElementType types[SMDSAbs_NbElementTypes] = { + SMDSAbs_Volume, SMDSAbs_Face, SMDSAbs_Edge, + SMDSAbs_0DElement, SMDSAbs_Ball, SMDSAbs_Node + }; + for ( int i = 0; i < SMDSAbs_NbElementTypes; ++i ) + if ( mesh->GetMeshInfo().NbElements( types[i] )) + { + type = types[i]; + break; + } + // put all elements in the vector + allElems.reserve( mesh->GetMeshInfo().NbElements( type )); + elemIt = mesh->elementsIterator( type ); + while ( elemIt->more() ) + allElems.push_back( elemIt->next()); + elemIt = elemSetIterator( allElems ); + } + else + { + type = (*theElements.begin())->GetType(); + elemIt = elemSetIterator( theElements ); + } + + // duplicate elements + + if ( type == SMDSAbs_Ball ) + { + SMDS_UnstructuredGrid* vtkGrid = mesh->getGrid(); + while ( elemIt->more() ) + { + const SMDS_MeshElement* elem = elemIt->next(); + if ( elem->GetType() != SMDSAbs_Ball ) + continue; + if (( elem = mesh->AddBall( elem->GetNode(0), + vtkGrid->GetBallDiameter( elem->getVtkId() )))) + myLastCreatedElems.Append( elem ); + } + } + else + { + vector< const SMDS_MeshNode* > nodes; + while ( elemIt->more() ) + { + const SMDS_MeshElement* elem = elemIt->next(); + if ( elem->GetType() != type ) + continue; + + nodes.assign( elem->begin_nodes(), elem->end_nodes() ); + + if ( type == SMDSAbs_Volume && elem->GetVtkType() == VTK_POLYHEDRON ) + { + std::vector quantities = + static_cast< const SMDS_VtkVolume* >( elem )->GetQuantities(); + elem = mesh->AddPolyhedralVolume( nodes, quantities ); + } + else + { + AddElement( nodes, type, elem->IsPoly() ); + elem = 0; // myLastCreatedElems is already filled + } + if ( elem ) + myLastCreatedElems.Append( elem ); + } + } +} + //================================================================================ /*! \brief Creates a hole in a mesh by doubling the nodes of some particular elements @@ -12144,7 +10969,8 @@ void SMESH_MeshEditor::CreateHoleSkin(double radius, double z = nodesCoords[i++]; gp_Pnt p = gp_Pnt(x, y ,z); gpnts.push_back(p); - MESSAGE("TopoDS_Vertex " << k++ << " " << p.X() << " " << p.Y() << " " << p.Z()); + MESSAGE("TopoDS_Vertex " << k << " " << p.X() << " " << p.Y() << " " << p.Z()); + k++; } } else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius @@ -12642,22 +11468,20 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, typedef vector TConnectivity; SMDS_ElemIteratorPtr eIt; - if (elements.empty()) - eIt = aMesh->elementsIterator(elemType); - else - eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() )); + if (elements.empty()) eIt = aMesh->elementsIterator(elemType); + else eIt = elemSetIterator( elements ); while (eIt->more()) { const SMDS_MeshElement* elem = eIt->next(); - const int iQuad = elem->IsQuadratic(); + const int iQuad = elem->IsQuadratic(); // ------------------------------------------------------------------------------------ // 1. For an elem, get present bnd elements and connectivities of missing bnd elements // ------------------------------------------------------------------------------------ vector presentBndElems; vector missingBndElems; - TConnectivity nodes; + TConnectivity nodes, elemNodes; if ( vTool.Set(elem, /*ignoreCentralNodes=*/true) ) // elem is a volume -------------- { vTool.SetExternalNormal(); @@ -12667,8 +11491,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, if ( !vTool.IsFreeFace(iface, &otherVol) && ( !aroundElements || elements.count( otherVol ))) continue; - const int nbFaceNodes = vTool.NbFaceNodes(iface); const SMDS_MeshNode** nn = vTool.GetFaceNodes(iface); + const int nbFaceNodes = vTool.NbFaceNodes (iface); if ( missType == SMDSAbs_Edge ) // boundary edges { nodes.resize( 2+iQuad ); @@ -12687,10 +11511,10 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, { nodes.clear(); for ( inode = 0; inode < nbFaceNodes; inode += 1+iQuad) - nodes.push_back( nn[inode] ); - if (iQuad) // add medium nodes + nodes.push_back( nn[inode] ); // add corner nodes + if (iQuad) for ( inode = 1; inode < nbFaceNodes; inode += 2) - nodes.push_back( nn[inode] ); + nodes.push_back( nn[inode] ); // add medium nodes int iCenter = vTool.GetCenterNodeIndex(iface); // for HEX27 if ( iCenter > 0 ) nodes.push_back( vTool.GetNodes()[ iCenter ] ); @@ -12718,22 +11542,24 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, } } } - else // elem is a face ------------------------------------------ + else if ( elem->GetType() == SMDSAbs_Face ) // elem is a face ------------------------ { avoidSet.clear(), avoidSet.insert( elem ); - int nbNodes = elem->NbCornerNodes(); - nodes.resize( 2 /*+ iQuad*/); - for ( int i = 0; i < nbNodes; i++ ) + elemNodes.assign( SMDS_MeshElement::iterator( elem->interlacedNodesElemIterator() ), + SMDS_MeshElement::iterator() ); + elemNodes.push_back( elemNodes[0] ); + nodes.resize( 2 + iQuad ); + const int nbLinks = elem->NbCornerNodes(); + for ( int i = 0, iN = 0; i < nbLinks; i++, iN += 1+iQuad ) { - nodes[0] = elem->GetNode(i); - nodes[1] = elem->GetNode((i+1)%nbNodes); - if ( FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet)) + nodes[0] = elemNodes[iN]; + nodes[1] = elemNodes[iN+1+iQuad]; + if ( SMESH_MeshAlgos::FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet)) continue; // not free link - //if ( iQuad ) - //nodes[2] = elem->GetNode( i + nbNodes ); + if ( iQuad ) nodes[2] = elemNodes[iN+1]; if ( const SMDS_MeshElement* edge = - aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/true)) + aMesh->FindElement(nodes,SMDSAbs_Edge,/*noMedium=*/false)) presentBndElems.push_back( edge ); else missingBndElems.push_back( nodes ); @@ -12837,10 +11663,8 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, // ----------------------- if ( toCopyElements && targetMesh != myMesh ) { - if (elements.empty()) - eIt = aMesh->elementsIterator(elemType); - else - eIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() )); + if (elements.empty()) eIt = aMesh->elementsIterator(elemType); + else eIt = elemSetIterator( elements ); while (eIt->more()) { const SMDS_MeshElement* elem = eIt->next();