X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=c62006a0eb2508d7ec9df9ba2f59518ba034dc4c;hp=f5bc661fbe1db0e0597e36aff91ac17280bab22c;hb=664ae5e0335821f9ca1ca863d3853841577ffd7c;hpb=85cd6043817a6218e3ff7c01b62f247d791a434f diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index f5bc661fb..c62006a0e 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 @@ -166,6 +162,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); @@ -527,12 +529,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]; @@ -545,53 +547,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]) { @@ -603,18 +594,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; } @@ -703,9 +694,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 @@ -718,32 +711,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; } @@ -765,30 +787,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; } //======================================================================= @@ -977,9 +997,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) @@ -1031,82 +1051,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; } @@ -1137,7 +1124,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 ); @@ -1186,8 +1173,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 ); @@ -1200,10 +1188,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]; @@ -1277,7 +1265,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; @@ -1297,13 +1286,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 ) { @@ -1317,70 +1305,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 @@ -1393,10 +1339,9 @@ 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; @@ -2377,18 +2322,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 @@ -2533,31 +2479,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]; @@ -2568,44 +2514,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]; @@ -2616,19 +2568,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] ); } } @@ -3587,39 +3546,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() ); } } } @@ -3646,7 +3600,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; } @@ -3901,7 +3855,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 @@ -4114,7 +4069,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++ ) @@ -4201,7 +4156,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; @@ -4227,9 +4182,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 @@ -4258,10 +4213,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 ); @@ -4277,6 +4237,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 @@ -4449,56 +4411,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 ) + // 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[ 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 ) - { - 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 } @@ -6184,1458 +6138,89 @@ void SMESH_MeshEditor::FindCoincidentNodes (TIDSortedNodeSet & theNodes, SMESH_OctreeNode::FindCoincidentNodes ( theNodes, &theGroupsOfNodes, theTolerance); } - //======================================================================= -/*! - * \brief Implementation of search for the node closest to point - */ +//function : SimplifyFace +//purpose : //======================================================================= -struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher +int SMESH_MeshEditor::SimplifyFace (const vector& faceNodes, + vector& poly_nodes, + vector& quantities) const { - //--------------------------------------------------------------------- - /*! - * \brief Constructor - */ - SMESH_NodeSearcherImpl( const SMESHDS_Mesh* theMesh ) - { - myMesh = ( SMESHDS_Mesh* ) theMesh; + int nbNodes = faceNodes.size(); - 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) ; + if (nbNodes < 3) + return 0; - // 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(); + set nodeSet; + + // get simple seq of nodes + //const SMDS_MeshNode* simpleNodes[ nbNodes ]; + vector simpleNodes( nbNodes ); + int iSimple = 0, nbUnique = 0; + + simpleNodes[iSimple++] = faceNodes[0]; + nbUnique++; + for (int iCur = 1; iCur < nbNodes; iCur++) { + if (faceNodes[iCur] != simpleNodes[iSimple - 1]) { + simpleNodes[iSimple++] = faceNodes[iCur]; + if (nodeSet.insert( faceNodes[iCur] ).second) + nbUnique++; } - 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() ); + int nbSimple = iSimple; + if (simpleNodes[nbSimple - 1] == simpleNodes[0]) { + nbSimple--; + iSimple--; } - //--------------------------------------------------------------------- - /*! - * \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() ); + if (nbUnique < 3) + return 0; + + // separate loops + int nbNew = 0; + bool foundLoop = (nbSimple > nbUnique); + while (foundLoop) { + foundLoop = false; + set loopSet; + for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) { + const SMDS_MeshNode* n = simpleNodes[iSimple]; + if (!loopSet.insert( n ).second) { + foundLoop = true; + + // separate loop + int iC = 0, curLast = iSimple; + for (; iC < curLast; iC++) { + if (simpleNodes[iC] == n) break; } - 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; + int loopLen = curLast - iC; + if (loopLen > 2) { + // create sub-element + nbNew++; + quantities.push_back(loopLen); + for (; iC < curLast; iC++) { + poly_nodes.push_back(simpleNodes[iC]); + } + } + // shift the rest nodes (place from the first loop position) + for (iC = curLast + 1; iC < nbSimple; iC++) { + simpleNodes[iC - loopLen] = simpleNodes[iC]; + } + nbSimple -= loopLen; + iSimple -= loopLen; } - } - return closestNode; - } - - //--------------------------------------------------------------------- - /*! - * \brief Destructor - */ - ~SMESH_NodeSearcherImpl() { delete myOctreeNode; } + } // for (iSimple = 0; iSimple < nbSimple; iSimple++) + } // while (foundLoop) - //--------------------------------------------------------------------- - /*! - * \brief Return the node tree - */ - const SMESH_OctreeNode* getTree() const { return myOctreeNode; } + if (iSimple > 2) { + nbNew++; + quantities.push_back(iSimple); + for (int i = 0; i < iSimple; i++) + poly_nodes.push_back(simpleNodes[i]); + } -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 ) - { - SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] ); - if ( !IsOut( &edge, point, tol )) - return false; - } - 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 ) - { - iClosest = i; - rClosest = r; - pClosest = pInt; - distClosest = intDist; - } - } - } - 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; - - // 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; - - bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0; - return covexCorner ? (out || out2) : (out && out2); - } - if ( element->GetType() == SMDSAbs_Edge ) // -------------------------------------------------- - { - // 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 ) - { - 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; -} - -//======================================================================= - -namespace -{ - // Position of a point relative to a segment - // . . - // . LEFT . - // . . - // VERTEX 1 o----ON-----> VERTEX 2 - // . . - // . RIGHT . - // . . - enum PositionName { POS_LEFT = 1, POS_VERTEX = 2, POS_RIGHT = 4, //POS_ON = 8, - POS_ALL = POS_LEFT | POS_RIGHT | POS_VERTEX }; - struct PointPos - { - PositionName _name; - int _index; // index of vertex or segment - - PointPos( PositionName n, int i=-1 ): _name(n), _index(i) {} - bool operator < (const PointPos& other ) const - { - if ( _name == other._name ) - return ( _index < 0 || other._index < 0 ) ? false : _index < other._index; - return _name < other._name; - } - }; - - //================================================================================ - /*! - * \brief Return of a point relative to a segment - * \param point2D - the point to analyze position of - * \param xyVec - end points of segments - * \param index0 - 0-based index of the first point of segment - * \param posToFindOut - flags of positions to detect - * \retval PointPos - point position - */ - //================================================================================ - - PointPos getPointPosition( const gp_XY& point2D, - const gp_XY* segEnds, - const int index0 = 0, - const int posToFindOut = POS_ALL) - { - const gp_XY& p1 = segEnds[ index0 ]; - const gp_XY& p2 = segEnds[ index0+1 ]; - const gp_XY grad = p2 - p1; - - if ( posToFindOut & POS_VERTEX ) - { - // check if the point2D is at "vertex 1" zone - gp_XY pp1[2] = { p1, gp_XY( p1.X() - grad.Y(), - p1.Y() + grad.X() ) }; - if ( getPointPosition( point2D, pp1, 0, POS_LEFT|POS_RIGHT )._name == POS_LEFT ) - return PointPos( POS_VERTEX, index0 ); - - // check if the point2D is at "vertex 2" zone - gp_XY pp2[2] = { p2, gp_XY( p2.X() - grad.Y(), - p2.Y() + grad.X() ) }; - if ( getPointPosition( point2D, pp2, 0, POS_LEFT|POS_RIGHT )._name == POS_RIGHT ) - return PointPos( POS_VERTEX, index0 + 1); - } - double edgeEquation = - ( point2D.X() - p1.X() ) * grad.Y() - ( point2D.Y() - p1.Y() ) * grad.X(); - return PointPos( edgeEquation < 0 ? POS_LEFT : POS_RIGHT, index0 ); - } -} - -//======================================================================= -/*! - * \brief Return minimal distance from a point to a face - * - * Currently we ignore non-planarity and 2nd order of face - */ -//======================================================================= - -double SMESH_MeshEditor::GetDistance( const SMDS_MeshFace* face, - const gp_Pnt& point ) -{ - double badDistance = -1; - if ( !face ) return badDistance; - - // coordinates of nodes (medium nodes, if any, ignored) - typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > TXyzIterator; - vector 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() ); - - // loop on segments of the face to analyze point position ralative to the face - set< PointPos > pntPosSet; - for ( size_t i = 1; i < xy.size(); ++i ) - { - PointPos pos = getPointPosition( point2D, &xy[0], i-1 ); - pntPosSet.insert( pos ); - } - - // compute distance - PointPos pos = *pntPosSet.begin(); - // cout << "Face " << face->GetID() << " DIST: "; - switch ( pos._name ) - { - case POS_LEFT: { - // point is most close to a segment - gp_Vec p0p1( point, xyz[ pos._index ] ); - gp_Vec p1p2( xyz[ pos._index ], xyz[ pos._index+1 ]); // segment vector - p1p2.Normalize(); - double projDist = p0p1 * p1p2; // distance projected to the segment - gp_Vec projVec = p1p2 * projDist; - gp_Vec distVec = p0p1 - projVec; - // cout << distVec.Magnitude() << ", SEG " << face->GetNode(pos._index)->GetID() - // << " - " << face->GetNodeWrap(pos._index+1)->GetID() << endl; - return distVec.Magnitude(); - } - case POS_RIGHT: { - // point is inside the face - double distToFacePlane = tmpPnt.Y(); - // cout << distToFacePlane << ", INSIDE " << endl; - return Abs( distToFacePlane ); - } - case POS_VERTEX: { - // point is most close to a node - gp_Vec distVec( point, xyz[ pos._index ]); - // cout << distVec.Magnitude() << " VERTEX " << face->GetNode(pos._index)->GetID() << endl; - return distVec.Magnitude(); - } - } - return badDistance; -} - -//======================================================================= -//function : SimplifyFace -//purpose : -//======================================================================= -int SMESH_MeshEditor::SimplifyFace (const vector faceNodes, - vector& poly_nodes, - vector& quantities) const -{ - int nbNodes = faceNodes.size(); - - if (nbNodes < 3) - return 0; - - set nodeSet; - - // get simple seq of nodes - //const SMDS_MeshNode* simpleNodes[ nbNodes ]; - vector simpleNodes( nbNodes ); - int iSimple = 0, nbUnique = 0; - - simpleNodes[iSimple++] = faceNodes[0]; - nbUnique++; - for (int iCur = 1; iCur < nbNodes; iCur++) { - if (faceNodes[iCur] != simpleNodes[iSimple - 1]) { - simpleNodes[iSimple++] = faceNodes[iCur]; - if (nodeSet.insert( faceNodes[iCur] ).second) - nbUnique++; - } - } - int nbSimple = iSimple; - if (simpleNodes[nbSimple - 1] == simpleNodes[0]) { - nbSimple--; - iSimple--; - } - - if (nbUnique < 3) - return 0; - - // separate loops - int nbNew = 0; - bool foundLoop = (nbSimple > nbUnique); - while (foundLoop) { - foundLoop = false; - set loopSet; - for (iSimple = 0; iSimple < nbSimple && !foundLoop; iSimple++) { - const SMDS_MeshNode* n = simpleNodes[iSimple]; - if (!loopSet.insert( n ).second) { - foundLoop = true; - - // separate loop - int iC = 0, curLast = iSimple; - for (; iC < curLast; iC++) { - if (simpleNodes[iC] == n) break; - } - int loopLen = curLast - iC; - if (loopLen > 2) { - // create sub-element - nbNew++; - quantities.push_back(loopLen); - for (; iC < curLast; iC++) { - poly_nodes.push_back(simpleNodes[iC]); - } - } - // shift the rest nodes (place from the first loop position) - for (iC = curLast + 1; iC < nbSimple; iC++) { - simpleNodes[iC - loopLen] = simpleNodes[iC]; - } - nbSimple -= loopLen; - iSimple -= loopLen; - } - } // for (iSimple = 0; iSimple < nbSimple; iSimple++) - } // while (foundLoop) - - if (iSimple > 2) { - nbNew++; - quantities.push_back(iSimple); - for (int i = 0; i < iSimple; i++) - poly_nodes.push_back(simpleNodes[i]); - } - - return nbNew; -} + return nbNew; +} //======================================================================= //function : MergeNodes @@ -8428,75 +7013,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 : @@ -8509,7 +7025,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 ); } //======================================================================= @@ -8625,7 +7141,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(); } @@ -9609,24 +8125,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(); @@ -9636,6 +8175,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 ) @@ -9657,7 +8202,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, break; default: NewElem = theHelper.AddPolygonalFace(nodes, id, theForce3d); - continue; } break; } @@ -9690,7 +8234,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; @@ -9710,6 +8254,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() ) { @@ -9725,19 +8270,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); @@ -9746,16 +8294,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(); @@ -9767,28 +8335,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 ) @@ -9809,6 +8391,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], @@ -9828,8 +8413,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); } } @@ -9867,36 +8453,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 ); } @@ -9934,27 +8520,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); @@ -9991,12 +8588,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); } } @@ -10155,7 +8759,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; @@ -10520,8 +9124,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() <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(); @@ -12681,8 +11287,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 ); @@ -12701,10 +11307,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 ] ); @@ -12732,22 +11338,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 );