X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=6088625bd52e036abcd0f6b5def5d80c6fb34c85;hb=949330b9879fbe09adda1a9919b7e375dee9d323;hp=e3fbc6602c1db1b06f17eabac1e6fbc201f43374;hpb=d4a710ce52f6e76786a7b3845e2f7975dc9a00b1;p=modules%2Fsmesh.git diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index e3fbc6602..6088625bd 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -20,7 +20,6 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -// SMESH SMESH : idl implementation based on 'SMESH' unit's classes // File : SMESH_MeshEditor.cxx // Created : Mon Apr 12 16:10:22 2004 // Author : Edward AGAPOV (eap) @@ -78,6 +77,7 @@ #include #include #include +#include #include #include #include @@ -97,6 +97,9 @@ #include #include +#include +#include + #define cast2Node(elem) static_cast( elem ) using namespace std; @@ -117,6 +120,19 @@ SMESH_MeshEditor::SMESH_MeshEditor( SMESH_Mesh* theMesh ) { } +//================================================================================ +/*! + * \brief Clears myLastCreatedNodes and myLastCreatedElems + */ +//================================================================================ + +void SMESH_MeshEditor::CrearLastCreated() +{ + myLastCreatedNodes.Clear(); + myLastCreatedElems.Clear(); +} + + //======================================================================= /*! * \brief Add element @@ -127,7 +143,8 @@ SMDS_MeshElement* SMESH_MeshEditor::AddElement(const vector & node, const SMDSAbs_ElementType type, const bool isPoly, - const int ID) + const int ID, + const double ballDiameter) { //MESSAGE("AddElement " < & node, else e = mesh->AddNode (node[0]->X(), node[0]->Y(), node[0]->Z()); break; + case SMDSAbs_Ball: + if ( ID >= 1 ) e = mesh->AddBallWithID(node[0], ballDiameter, ID); + else e = mesh->AddBall (node[0], ballDiameter); + break; + default:; } if ( e ) myLastCreatedElems.Append( e ); @@ -380,6 +402,44 @@ int SMESH_MeshEditor::Remove (const list< int >& theIDs, return removed; } +//================================================================================ +/*! + * \brief Create 0D elements on all nodes of the given object except those + * nodes on which a 0D element already exists. + * \param elements - Elements on whose nodes to create 0D elements; if empty, + * the all mesh is treated + * \param all0DElems - returns all 0D elements found or created on nodes of \a elements + */ +//================================================================================ + +void SMESH_MeshEditor::Create0DElementsOnAllNodes( const TIDSortedElemSet& elements, + TIDSortedElemSet& all0DElems ) +{ + typedef SMDS_SetIterator TSetIterator; + SMDS_ElemIteratorPtr elemIt; + if ( elements.empty() ) + elemIt = GetMeshDS()->elementsIterator( SMDSAbs_Node ); + else + elemIt = SMDS_ElemIteratorPtr( new TSetIterator( elements.begin(), elements.end() )); + + while ( elemIt->more() ) + { + const SMDS_MeshElement* e = elemIt->next(); + SMDS_ElemIteratorPtr nodeIt = e->nodesIterator(); + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = cast2Node( nodeIt->next() ); + SMDS_ElemIteratorPtr it0D = n->GetInverseElementIterator( SMDSAbs_0DElement ); + if ( it0D->more() ) + all0DElems.insert( it0D->next() ); + else { + myLastCreatedElems.Append( GetMeshDS()->Add0DElement( n )); + all0DElems.insert( myLastCreatedElems.Last() ); + } + } + } +} + //======================================================================= //function : FindShape //purpose : Return an index of the shape theElem is on @@ -487,10 +547,10 @@ static void ShiftNodesQuadTria(const SMDS_MeshNode* aNodes[]) //======================================================================= //function : edgeConnectivity -//purpose : auxilary +//purpose : auxilary // return number of the edges connected with the theNode. // if theEdges has connections with the other type of the -// elements, return -1 +// elements, return -1 //======================================================================= static int nbEdgeConnectivity(const SMDS_MeshNode* theNode) { @@ -1051,6 +1111,97 @@ bool SMESH_MeshEditor::Reorient (const SMDS_MeshElement * theElem) return false; } +//================================================================================ +/*! + * \brief Reorient faces. + * \param theFaces - the faces to reorient. If empty the whole mesh is meant + * \param theDirection - desired direction of normal of \a theFace + * \param theFace - one of \a theFaces that sould be orientated according to + * \a theDirection and whose orientation defines orientation of other faces + * \return number of reoriented faces. + */ +//================================================================================ + +int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, + const gp_Dir& theDirection, + const SMDS_MeshElement * theFace) +{ + int nbReori = 0; + if ( !theFace || theFace->GetType() != SMDSAbs_Face ) return nbReori; + + if ( theFaces.empty() ) + { + SMDS_FaceIteratorPtr fIt = GetMeshDS()->facesIterator(/*idInceasingOrder=*/true); + while ( fIt->more() ) + theFaces.insert( theFaces.end(), fIt->next() ); + } + + // orient theFace according to theDirection + gp_XYZ normal; + SMESH_Algo::FaceNormal( theFace, normal, /*normalized=*/false ); + if ( normal * theDirection.XYZ() < 0 ) + nbReori += Reorient( theFace ); + + // Orient other faces + + set< const SMDS_MeshElement* > startFaces; + TIDSortedElemSet avoidSet; + set< SMESH_TLink > checkedLinks; + pair< set< SMESH_TLink >::iterator, bool > linkIt_isNew; + + if ( theFaces.size() > 1 )// leave 1 face to prevent finding not selected faces + theFaces.erase( theFace ); + startFaces.insert( theFace ); + + set< const SMDS_MeshElement* >::iterator startFace = startFaces.begin(); + while ( startFace != startFaces.end() ) + { + theFace = *startFace; + const int nbNodes = theFace->NbCornerNodes(); + + avoidSet.clear(); + avoidSet.insert(theFace); + + NLink link( theFace->GetNode( 0 ), 0 ); + for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace + { + link.second = theFace->GetNode(( i+1 ) % nbNodes ); + linkIt_isNew = checkedLinks.insert( link ); + if ( !linkIt_isNew.second ) + { + // link has already been checked and won't be encountered more + // if the group (theFaces) is manifold + checkedLinks.erase( linkIt_isNew.first ); + } + else + { + int nodeInd1, nodeInd2; + const SMDS_MeshElement* otherFace = FindFaceInSet( link.first, link.second, + theFaces, avoidSet, + & nodeInd1, & nodeInd2); + if ( otherFace && otherFace != theFace) + { + // link must be reversed in otherFace if orientation ot otherFace + // is same as that of theFace + if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 ) + { + // cout << "Reorient " << otherFace->GetID() << " near theFace=" <GetID() + // << " \tlink( " << link.first->GetID() << " " << link.second->GetID() << endl; + nbReori += Reorient( otherFace ); + } + startFaces.insert( otherFace ); + if ( theFaces.size() > 1 ) // leave 1 face to prevent finding not selected faces + theFaces.erase( otherFace ); + } + } + std::swap( link.first, link.second ); + } + startFaces.erase( startFace ); + startFace = startFaces.begin(); + } + return nbReori; +} + //======================================================================= //function : getBadRate //purpose : @@ -1116,7 +1267,8 @@ bool SMESH_MeshEditor::QuadToTri (TIDSortedElemSet & theElems, 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 ) { // tr1 + tr2 is better newElem1 = aMesh->AddFace( aNodes[2], aNodes[3], aNodes[0] ); @@ -1250,7 +1402,8 @@ int SMESH_MeshEditor::BestSplit (const SMDS_MeshElement* theQuad, SMDS_FaceOfNodes tr3 ( aNodes[1], aNodes[2], aNodes[3] ); SMDS_FaceOfNodes tr4 ( aNodes[3], aNodes[0], aNodes[1] ); aBadRate2 = getBadRate( &tr3, theCrit ) + getBadRate( &tr4, theCrit ); - + // 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) // tr1 + tr2 is better return 1; // diagonal 1-3 @@ -1641,7 +1794,7 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1); SMESHDS_SubMesh* fSubMesh = 0;//subMesh; - + SMESH_SequenceOfElemPtr newNodes, newElems; // map face of volume to it's baricenrtic node @@ -2754,7 +2907,7 @@ void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode, const SMDS_MeshElement* elem = elemIt->next(); if(elem->GetType() == SMDSAbs_0DElement) continue; - + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); if ( elem->GetType() == SMDSAbs_Volume ) { @@ -2962,7 +3115,7 @@ void SMESH_MeshEditor::Smooth (TIDSortedElemSet & theElems, SMDS_FaceIteratorPtr fIt = aMesh->facesIterator(); while ( fIt->more() ) { const SMDS_MeshElement* face = fIt->next(); - theElems.insert( face ); + theElems.insert( theElems.end(), face ); } } // get all face ids theElems are on @@ -3481,9 +3634,9 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem, //MESSAGE("sweepElement " << nbSteps); SMESHDS_Mesh* aMesh = GetMeshDS(); - const int nbNodes = elem->NbNodes(); + const int nbNodes = elem->NbNodes(); const int nbCorners = elem->NbCornerNodes(); - SMDSAbs_EntityType baseType = elem->GetEntityType(); /* it can change in case of + SMDSAbs_EntityType baseType = elem->GetEntityType(); /* it can change in case of polyhedron creation !!! */ // Loop on elem nodes: // find new nodes and detect same nodes indices @@ -3761,7 +3914,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem, midlNod[0], midlNod[1], midlNod[2], midlNod[3]); } else if(nbSame==1) { - // ---> pyramid + pentahedron - can not be created since it is needed + // ---> pyramid + pentahedron - can not be created since it is needed // additional middle node at the center of face INFOS( " Sweep for face " << elem->GetID() << " can not be created" ); return; @@ -3827,6 +3980,9 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem, } break; } + case SMDSEntity_Ball: + return; + default: break; } @@ -3928,6 +4084,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, { const SMDS_MeshNode* node = static_cast( nList->first ); + if ( newElemsMap.count( node )) + continue; // node was extruded into edge SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator(); int nbInitElems = 0; const SMDS_MeshElement* el = 0; @@ -4719,7 +4877,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, list< list > LLPPs; int startNid = theN1->GetID(); TColStd_MapOfInteger UsedNums; - + int NbEdges = Edges.Length(); int i = 1; for(; i<=NbEdges; i++) { @@ -4865,13 +5023,12 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, const SMDS_MeshElement* currentElem = NULL; int totalNbEdges = theTrack->NbEdges(); SMDS_ElemIteratorPtr nIt; - bool isClosed = false; //check start node if( !theTrack->GetMeshDS()->Contains(theN1) ) { return EXTR_BAD_STARTING_NODE; } - + conn = nbEdgeConnectivity(theN1); if(conn > 2) return EXTR_PATH_NOT_EDGE; @@ -4886,24 +5043,23 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, if(currentNode == prevNode) currentNode = static_cast(nIt->next()); aNodesList.push_back(currentNode); - } else { + } else { nIt = currentElem->nodesIterator(); while( nIt->more() ) { currentNode = static_cast(nIt->next()); if(currentNode == prevNode) currentNode = static_cast(nIt->next()); aNodesList.push_back(currentNode); - + //case of the closed mesh if(currentNode == theN1) { nbEdges++; - isClosed = true; break; } conn = nbEdgeConnectivity(currentNode); if(conn > 2) { - return EXTR_PATH_NOT_EDGE; + return EXTR_PATH_NOT_EDGE; }else if( conn == 1 && nbEdges > 0 ) { //End of the path nbEdges++; @@ -4919,8 +5075,8 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, nbEdges++; } } - } - + } + if(nbEdges != totalNbEdges) return EXTR_PATH_NOT_EDGE; @@ -4932,7 +5088,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, x1 = aNodesList[i-1]->X();x2 = aNodesList[i]->X(); y1 = aNodesList[i-1]->Y();y2 = aNodesList[i]->Y(); z1 = aNodesList[i-1]->Z();z2 = aNodesList[i]->Z(); - TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2)); + TopoDS_Edge e = BRepBuilderAPI_MakeEdge(gp_Pnt(x1,y1,z1),gp_Pnt(x2,y2,z2)); list LPP; aPrms.clear(); MakeEdgePathPoints(aPrms, e, (aNodesList[i-1]->GetID()==startNid), LPP); @@ -4971,7 +5127,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, fullList.pop_back(); } fullList.push_back(PP1); - + } // Sub-shape for the Pattern must be an Edge or Wire else if( aS.ShapeType() == TopAbs_EDGE ) { aTrackEdge = TopoDS::Edge( aS ); @@ -5602,139 +5758,154 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems, for ( itElem = theElems.begin(); itElem != theElems.end(); itElem++ ) { const SMDS_MeshElement* elem = *itElem; - if ( !elem || elem->GetType() == SMDSAbs_Node ) - continue; + if ( !elem ) continue; - int nbNodes = elem->NbNodes(); - int elemType = elem->GetType(); + SMDSAbs_GeometryType geomType = elem->GetGeomType(); + int nbNodes = elem->NbNodes(); + if ( geomType == SMDSGeom_NONE ) continue; // node - if (elem->IsPoly()) { + switch ( geomType ) { - // polygon or polyhedral volume - switch ( elemType ) { - case SMDSAbs_Face: - { - vector poly_nodes (nbNodes); - int iNode = 0; - SMDS_ElemIteratorPtr itN = elem->nodesIterator(); - while (itN->more()) { - const SMDS_MeshNode* node = - static_cast(itN->next()); - TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); - if (nodeMapIt == nodeMap.end()) - break; // not all nodes transformed - if (needReverse) { - // reverse mirrored faces and volumes - poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second; - } else { - poly_nodes[iNode] = (*nodeMapIt).second; - } - iNode++; + case SMDSGeom_POLYGON: // ---------------------- polygon + { + vector poly_nodes (nbNodes); + int iNode = 0; + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while (itN->more()) { + const SMDS_MeshNode* node = + static_cast(itN->next()); + TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); + if (nodeMapIt == nodeMap.end()) + break; // not all nodes transformed + if (needReverse) { + // reverse mirrored faces and volumes + poly_nodes[nbNodes - iNode - 1] = (*nodeMapIt).second; + } else { + poly_nodes[iNode] = (*nodeMapIt).second; } - if ( iNode != nbNodes ) - continue; // not all nodes transformed + iNode++; + } + if ( iNode != nbNodes ) + continue; // not all nodes transformed - if ( theTargetMesh ) { - myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes)); - srcElems.Append( elem ); - } - else if ( theCopy ) { - myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes)); - srcElems.Append( elem ); - } - else { - aMesh->ChangePolygonNodes(elem, poly_nodes); - } + if ( theTargetMesh ) { + myLastCreatedElems.Append(aTgtMesh->AddPolygonalFace(poly_nodes)); + srcElems.Append( elem ); } - break; - case SMDSAbs_Volume: - { - const SMDS_VtkVolume* aPolyedre = - dynamic_cast( elem ); - if (!aPolyedre) { - MESSAGE("Warning: bad volumic element"); - continue; - } + else if ( theCopy ) { + myLastCreatedElems.Append(aMesh->AddPolygonalFace(poly_nodes)); + srcElems.Append( elem ); + } + else { + aMesh->ChangePolygonNodes(elem, poly_nodes); + } + } + break; - vector poly_nodes; poly_nodes.reserve( nbNodes ); - vector quantities; + case SMDSGeom_POLYHEDRA: // ------------------ polyhedral volume + { + const SMDS_VtkVolume* aPolyedre = + dynamic_cast( elem ); + if (!aPolyedre) { + MESSAGE("Warning: bad volumic element"); + continue; + } - bool allTransformed = true; - int nbFaces = aPolyedre->NbFaces(); - for (int iface = 1; iface <= nbFaces && allTransformed; iface++) { - int nbFaceNodes = aPolyedre->NbFaceNodes(iface); - for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) { - const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode); - TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); - if (nodeMapIt == nodeMap.end()) { - allTransformed = false; // not all nodes transformed - } else { - poly_nodes.push_back((*nodeMapIt).second); - } - if ( needReverse && allTransformed ) - std::reverse( poly_nodes.end() - nbFaceNodes, poly_nodes.end() ); + vector poly_nodes; poly_nodes.reserve( nbNodes ); + vector quantities; quantities.reserve( nbNodes ); + + bool allTransformed = true; + int nbFaces = aPolyedre->NbFaces(); + for (int iface = 1; iface <= nbFaces && allTransformed; iface++) { + int nbFaceNodes = aPolyedre->NbFaceNodes(iface); + for (int inode = 1; inode <= nbFaceNodes && allTransformed; inode++) { + const SMDS_MeshNode* node = aPolyedre->GetFaceNode(iface, inode); + TNodeNodeMap::iterator nodeMapIt = nodeMap.find(node); + if (nodeMapIt == nodeMap.end()) { + allTransformed = false; // not all nodes transformed + } else { + poly_nodes.push_back((*nodeMapIt).second); } - quantities.push_back(nbFaceNodes); + if ( needReverse && allTransformed ) + std::reverse( poly_nodes.end() - nbFaceNodes, poly_nodes.end() ); } - if ( !allTransformed ) - continue; // not all nodes transformed + quantities.push_back(nbFaceNodes); + } + if ( !allTransformed ) + continue; // not all nodes transformed - if ( theTargetMesh ) { - myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities)); - srcElems.Append( elem ); - } - else if ( theCopy ) { - myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities)); - srcElems.Append( elem ); - } - else { - aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities); - } + if ( theTargetMesh ) { + myLastCreatedElems.Append(aTgtMesh->AddPolyhedralVolume(poly_nodes, quantities)); + srcElems.Append( elem ); + } + else if ( theCopy ) { + myLastCreatedElems.Append(aMesh->AddPolyhedralVolume(poly_nodes, quantities)); + srcElems.Append( elem ); + } + else { + aMesh->ChangePolyhedronNodes(elem, poly_nodes, quantities); } - break; - default:; } - continue; + break; - } // elem->isPoly() + case SMDSGeom_BALL: // -------------------- Ball + { + if ( !theCopy && !theTargetMesh ) continue; - // Regular elements + TNodeNodeMap::iterator nodeMapIt = nodeMap.find( elem->GetNode(0) ); + if (nodeMapIt == nodeMap.end()) + continue; // not all nodes transformed - while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() ); - const std::vector& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType() ); - const std::vector& i = needReverse ? iRev : iForw; + double diameter = static_cast(elem)->GetDiameter(); + if ( theTargetMesh ) { + myLastCreatedElems.Append(aTgtMesh->AddBall( nodeMapIt->second, diameter )); + srcElems.Append( elem ); + } + else { + myLastCreatedElems.Append(aMesh->AddBall( nodeMapIt->second, diameter )); + srcElems.Append( elem ); + } + } + break; - // find transformed nodes - vector nodes(nbNodes); - int iNode = 0; - SMDS_ElemIteratorPtr itN = elem->nodesIterator(); - while ( itN->more() ) { - const SMDS_MeshNode* node = - static_cast( itN->next() ); - TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node ); - if ( nodeMapIt == nodeMap.end() ) - break; // not all nodes transformed - nodes[ i [ iNode++ ]] = (*nodeMapIt).second; - } - if ( iNode != nbNodes ) - continue; // not all nodes transformed + default: // ----------------------- Regular elements - if ( theTargetMesh ) { - if ( SMDS_MeshElement* copy = - targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) { - myLastCreatedElems.Append( copy ); - srcElems.Append( elem ); + while ( iForw.size() < nbNodes ) iForw.push_back( iForw.size() ); + const std::vector& iRev = SMDS_MeshCell::reverseSmdsOrder( elem->GetEntityType() ); + const std::vector& i = needReverse ? iRev : iForw; + + // find transformed nodes + vector nodes(nbNodes); + int iNode = 0; + SMDS_ElemIteratorPtr itN = elem->nodesIterator(); + while ( itN->more() ) { + const SMDS_MeshNode* node = + static_cast( itN->next() ); + TNodeNodeMap::iterator nodeMapIt = nodeMap.find( node ); + if ( nodeMapIt == nodeMap.end() ) + break; // not all nodes transformed + nodes[ i [ iNode++ ]] = (*nodeMapIt).second; } - } - else if ( theCopy ) { - if ( AddElement( nodes, elem->GetType(), elem->IsPoly() )) - srcElems.Append( elem ); - } - else { - // reverse element as it was reversed by transformation - if ( nbNodes > 2 ) - aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes ); - } + if ( iNode != nbNodes ) + continue; // not all nodes transformed + + if ( theTargetMesh ) { + if ( SMDS_MeshElement* copy = + targetMeshEditor.AddElement( nodes, elem->GetType(), elem->IsPoly() )) { + myLastCreatedElems.Append( copy ); + srcElems.Append( elem ); + } + } + else if ( theCopy ) { + if ( AddElement( nodes, elem->GetType(), elem->IsPoly() )) + srcElems.Append( elem ); + } + else { + // reverse element as it was reversed by transformation + if ( nbNodes > 2 ) + aMesh->ChangeElementNodes( elem, &nodes[0], nbNodes ); + } + } // switch ( geomType ) } // loop on elements @@ -5768,26 +5939,32 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, // Sort existing groups by types and collect their names // to store an old group and a generated new one - typedef pair< SMESHDS_GroupBase*, SMDS_MeshGroup* > TOldNewGroup; + typedef pair< SMESHDS_GroupBase*, SMESHDS_Group* > TOldNewGroup; vector< list< TOldNewGroup > > groupsByType( SMDSAbs_NbElementTypes ); + vector< TOldNewGroup* > orderedOldNewGroups; // in order of old groups // group names set< string > groupNames; - // - SMDS_MeshGroup* nullNewGroup = (SMDS_MeshGroup*) 0; + SMESH_Mesh::GroupIteratorPtr groupIt = GetMesh()->GetGroups(); - while ( groupIt->more() ) { + if ( !groupIt->more() ) return newGroupIDs; + + int newGroupID = mesh->GetGroupIds().back()+1; + while ( groupIt->more() ) + { SMESH_Group * group = groupIt->next(); if ( !group ) continue; SMESHDS_GroupBase* groupDS = group->GetGroupDS(); if ( !groupDS || groupDS->IsEmpty() ) continue; groupNames.insert( group->GetName() ); groupDS->SetStoreName( group->GetName() ); - groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, nullNewGroup )); + SMESHDS_Group* newGroup = new SMESHDS_Group( newGroupID++, mesh->GetMeshDS(), + groupDS->GetType() ); + groupsByType[ groupDS->GetType() ].push_back( make_pair( groupDS, newGroup )); + orderedOldNewGroups.push_back( & groupsByType[ groupDS->GetType() ].back() ); } - // Groups creation + // Loop on nodes and elements to add them in new groups - // loop on nodes and elements for ( int isNodes = 0; isNodes < 2; ++isNodes ) { const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens; @@ -5804,7 +5981,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, continue; } list< TOldNewGroup > & groupsOldNew = groupsByType[ sourceElem->GetType() ]; - if ( groupsOldNew.empty() ) { + if ( groupsOldNew.empty() ) { // no groups of this type at all while ( iElem < gens.Length() && gens( iElem+1 ) == sourceElem ) ++iElem; // skip all elements made by sourceElem continue; @@ -5818,58 +5995,56 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, if ( const SMDS_MeshElement* resElem = elems( ++iElem )) if ( resElem != sourceElem ) resultElems.push_back( resElem ); - // do not generate element groups from node ones -// if ( sourceElem->GetType() == SMDSAbs_Node && -// elems( iElem )->GetType() != SMDSAbs_Node ) -// continue; // add resultElems to groups made by ones the sourceElem belongs to list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end(); for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew ) { SMESHDS_GroupBase* oldGroup = gOldNew->first; - if ( oldGroup->Contains( sourceElem )) // sourceElem in oldGroup + if ( oldGroup->Contains( sourceElem )) // sourceElem is in oldGroup { - SMDS_MeshGroup* & newGroup = gOldNew->second; - if ( !newGroup )// create a new group - { - // make a name - string name = oldGroup->GetStoreName(); - if ( !targetMesh ) { - name += "_"; - name += postfix; - int nb = 0; - while ( !groupNames.insert( name ).second ) // name exists - { - if ( nb == 0 ) { - name += "_1"; - } - else { - TCollection_AsciiString nbStr(nb+1); - name.resize( name.rfind('_')+1 ); - name += nbStr.ToCString(); - } - ++nb; - } - } - // make a group - int id; - SMESH_Group* group = mesh->AddGroup( resultElems.back()->GetType(), - name.c_str(), id ); - SMESHDS_Group* groupDS = static_cast(group->GetGroupDS()); - newGroup = & groupDS->SMDSGroup(); - newGroupIDs->push_back( id ); - } - // fill in a new group + SMDS_MeshGroup & newGroup = gOldNew->second->SMDSGroup(); list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt; for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt ) - newGroup->Add( *resElemIt ); + newGroup.Add( *resElemIt ); } } } // loop on created elements }// loop on nodes and elements + // Create new SMESH_Groups from SMESHDS_Groups and remove empty SMESHDS_Groups + + for ( size_t i = 0; i < orderedOldNewGroups.size(); ++i ) + { + SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->first; + SMESHDS_Group* newGroupDS = orderedOldNewGroups[i]->second; + if ( newGroupDS->IsEmpty() ) + { + mesh->GetMeshDS()->RemoveGroup( newGroupDS ); + } + else + { + // make a name + string name = oldGroupDS->GetStoreName(); + if ( !targetMesh ) { + name += "_"; + name += postfix; + int nb = 1; + while ( !groupNames.insert( name ).second ) // name exists + name = SMESH_Comment( oldGroupDS->GetStoreName() ) << "_" << postfix << "_" << nb++; + } + newGroupDS->SetStoreName( name.c_str() ); + + // make a SMESH_Groups + mesh->AddGroup( newGroupDS ); + newGroupIDs->push_back( newGroupDS->GetID() ); + + // set group type + newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() ); + } + } + return newGroupIDs; } @@ -5981,7 +6156,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher } else if ( tree->NbNodes() ) // put a tree to the treeMap { - const Bnd_B3d& box = tree->getBox(); + 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 @@ -5993,7 +6168,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher TDistTreeMap::iterator sqDist_tree = treeMap.begin(); if ( treeMap.size() > 5 ) { SMESH_OctreeNode* closestTree = sqDist_tree->second; - const Bnd_B3d& box = closestTree->getBox(); + const Bnd_B3d& box = *closestTree->getBox(); double limit = sqrt( sqDist_tree->first ) + sqrt ( box.SquareExtent() ); sqLimit = limit * limit; } @@ -6043,7 +6218,7 @@ private: */ //======================================================================= -SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher() +SMESH_NodeSearcher* SMESH_MeshEditor::GetNodeSearcher() { return new SMESH_NodeSearcherImpl( GetMeshDS() ); } @@ -6065,16 +6240,22 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() { public: - ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), double tolerance = NodeRadius ); - void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems); + ElementBndBoxTree(const SMDS_Mesh& mesh, + SMDSAbs_ElementType elemType, + SMDS_ElemIteratorPtr theElemIt = SMDS_ElemIteratorPtr(), + double tolerance = NodeRadius ); + void getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems ); void getElementsNearLine ( const gp_Ax1& line, TIDSortedElemSet& foundElems); + void getElementsInSphere ( const gp_XYZ& center, + const double radius, TIDSortedElemSet& foundElems); + size_t getSize() { return std::max( _size, _elements.size() ); } ~ElementBndBoxTree(); protected: - ElementBndBoxTree() {} - SMESH_Octree* allocateOctreeChild() const { return new ElementBndBoxTree; } - void buildChildrenData(); - Bnd_B3d* buildRootBox(); + 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 @@ -6084,6 +6265,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() ElementBox(const SMDS_MeshElement* elem, double tolerance); }; vector< ElementBox* > _elements; + size_t _size; }; //================================================================================ @@ -6093,7 +6275,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() //================================================================================ ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance) - :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. )) + :SMESH_Octree( new SMESH_TreeLimit( MaxLevel, /*minSize=*/0. )) { int nbElems = mesh.GetMeshInfo().NbElements( elemType ); _elements.reserve( nbElems ); @@ -6102,10 +6284,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() while ( elemIt->more() ) _elements.push_back( new ElementBox( elemIt->next(),tolerance )); - if ( _elements.size() > MaxNbElemsInLeaf ) - compute(); - else - myIsLeaf = true; + compute(); } //================================================================================ @@ -6147,7 +6326,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() { for (int j = 0; j < 8; j++) { - if ( !_elements[i]->IsOut( myChildren[j]->getBox() )) + if ( !_elements[i]->IsOut( *myChildren[j]->getBox() )) { _elements[i]->_refCount++; ((ElementBndBoxTree*)myChildren[j])->_elements.push_back( _elements[i]); @@ -6155,7 +6334,8 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() } _elements[i]->_refCount--; } - _elements.clear(); + _size = _elements.size(); + SMESHUtils::FreeVector( _elements ); // = _elements.clear() + free memory for (int j = 0; j < 8; j++) { @@ -6164,7 +6344,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() child->myIsLeaf = true; if ( child->_elements.capacity() - child->_elements.size() > 1000 ) - child->_elements.resize( child->_elements.size() ); // compact + SMESHUtils::CompactVector( child->_elements ); } } @@ -6177,7 +6357,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() void ElementBndBoxTree::getElementsNearPoint( const gp_Pnt& point, TIDSortedElemSet& foundElems) { - if ( level() && getBox().IsOut( point.XYZ() )) + if ( getBox()->IsOut( point.XYZ() )) return; if ( isLeaf() ) @@ -6202,7 +6382,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() void ElementBndBoxTree::getElementsNearLine( const gp_Ax1& line, TIDSortedElemSet& foundElems) { - if ( level() && getBox().IsOut( line )) + if ( getBox()->IsOut( line )) return; if ( isLeaf() ) @@ -6218,6 +6398,32 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() } } + //================================================================================ + /*! + * \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 @@ -6230,7 +6436,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() _refCount = 1; SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); while ( nIt->more() ) - Add( SMESH_TNodeXYZ( cast2Node( nIt->next() ))); + Add( SMESH_TNodeXYZ( nIt->next() )); Enlarge( tolerance ); } @@ -6265,6 +6471,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher 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, @@ -6375,12 +6583,12 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& lin 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) )); + SMESH_TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); anExtCC.Init( lineCurve, edge); if ( anExtCC.NbExtrema() > 0 && anExtCC.LowerDistance() <= tol) { @@ -6429,7 +6637,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF while (( f = SMESH_MeshEditor::FindFaceInSet(link.node1(), link.node2(), emptySet, faces ))) faces.insert( f ); - // select another outer face among the found + // select another outer face among the found const SMDS_MeshElement* outerFace2 = 0; if ( faces.size() == 2 ) { @@ -6489,7 +6697,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF // 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 { @@ -6502,7 +6710,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF * \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 and 0D elements + * 'ALL' type means elements of any type excluding nodes, balls and 0D elements */ //======================================================================= @@ -6516,7 +6724,7 @@ FindElementsByPoint(const gp_Pnt& point, double tolerance = getTolerance(); // ================================================================================= - if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement ) + if ( type == SMDSAbs_Node || type == SMDSAbs_0DElement || type == SMDSAbs_Ball) { if ( !_nodeSearcher ) _nodeSearcher = new SMESH_NodeSearcherImpl( _mesh ); @@ -6533,7 +6741,7 @@ FindElementsByPoint(const gp_Pnt& point, } else { - SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( SMDSAbs_0DElement ); + SMDS_ElemIteratorPtr elemIt = closeNode->GetInverseElementIterator( type ); while ( elemIt->more() ) foundElements.push_back( elemIt->next() ); } @@ -6550,12 +6758,103 @@ FindElementsByPoint(const gp_Pnt& point, _ebbTree->getElementsNearPoint( point, suspectElems ); TIDSortedElemSet::iterator elem = suspectElems.begin(); for ( ; elem != suspectElems.end(); ++elem ) - if ( !SMESH_MeshEditor::isOut( *elem, point, tolerance )) + if ( !SMESH_MeshEditor::IsOut( *elem, point, tolerance )) foundElements.push_back( *elem ); } return foundElements.size(); } +//======================================================================= +/*! + * \brief Find an element of given type most close to the given point + * + * WARNING: Only face search is implemeneted so far + */ +//======================================================================= + +const SMDS_MeshElement* +SMESH_ElementSearcherImpl::FindClosestTo( const gp_Pnt& point, + SMDSAbs_ElementType type ) +{ + const SMDS_MeshElement* closestElem = 0; + + if ( type == SMDSAbs_Face ) + { + if ( !_ebbTree || _elementType != type ) + { + if ( _ebbTree ) delete _ebbTree; + _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt ); + } + TIDSortedElemSet suspectElems; + _ebbTree->getElementsNearPoint( point, suspectElems ); + + if ( suspectElems.empty() && _ebbTree->maxSize() > 0 ) + { + gp_Pnt boxCenter = 0.5 * ( _ebbTree->getBox()->CornerMin() + + _ebbTree->getBox()->CornerMax() ); + double radius; + if ( _ebbTree->getBox()->IsOut( point.XYZ() )) + radius = point.Distance( boxCenter ) - 0.5 * _ebbTree->maxSize(); + else + radius = _ebbTree->maxSize() / pow( 2., _ebbTree->getHeight()) / 2; + while ( suspectElems.empty() ) + { + _ebbTree->getElementsInSphere( point.XYZ(), radius, suspectElems ); + radius *= 1.1; + } + } + double minDist = std::numeric_limits::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 @@ -6609,7 +6908,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) else if ( ! intersection.IsParallel() && intersection.NbPoints() > 0 ) { gp_Pnt intersectionPoint = intersection.Point(1); - if ( !SMESH_MeshEditor::isOut( *face, intersectionPoint, tolerance )) + if ( !SMESH_MeshEditor::IsOut( *face, intersectionPoint, tolerance )) u2inters.insert(make_pair( intersection.ParamOnConic(1), TInters( *face, fNorm ))); } } @@ -6617,7 +6916,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) int nbInter = u2inters.size(); if ( nbInter == 0 ) - return TopAbs_OUT; + return TopAbs_OUT; double f = u2inters.begin()->first, l = u2inters.rbegin()->first; if ( nbInter == 1 ) // not closed mesh @@ -6751,7 +7050,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) return TopAbs_ON; if ( nbIntBeforePoint == 0 || nbIntAfterPoint == 0) - return TopAbs_OUT; + return TopAbs_OUT; if ( nbIntBeforePoint + nbIntAfterPoint == 1 ) // not closed mesh return fabs( f ) < tolerance ? TopAbs_ON : TopAbs_UNKNOWN; @@ -6827,7 +7126,7 @@ SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher(SMDS_ElemIteratorPtr */ //======================================================================= -bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol ) +bool SMESH_MeshEditor::IsOut( const SMDS_MeshElement* element, const gp_Pnt& point, double tol ) { if ( element->GetType() == SMDSAbs_Volume) { @@ -6874,7 +7173,7 @@ bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& poi for ( i = 0; i < nbNodes; ++i ) { SMDS_LinearEdge edge( nodeList[i], nodeList[i+1] ); - if ( !isOut( &edge, point, tol )) + if ( !IsOut( &edge, point, tol )) return false; } return true; @@ -6981,22 +7280,186 @@ bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& poi } //======================================================================= -//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; +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 - set nodeSet; + 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; + } + }; - // get simple seq of nodes - //const SMDS_MeshNode* simpleNodes[ nbNodes ]; + //================================================================================ + /*! + * \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; @@ -7746,32 +8209,29 @@ private: //purpose : Return list of group of elements built on the same nodes. // Search among theElements or in the whole mesh if theElements is empty //======================================================================= -void SMESH_MeshEditor::FindEqualElements(set & theElements, - TListOfListOfElementsID & theGroupsOfElementsID) + +void SMESH_MeshEditor::FindEqualElements(TIDSortedElemSet & theElements, + TListOfListOfElementsID & theGroupsOfElementsID) { myLastCreatedElems.Clear(); myLastCreatedNodes.Clear(); - typedef set TElemsSet; typedef map< SortableElement, int > TMapOfNodeSet; typedef list TGroupOfElems; - TElemsSet elems; if ( theElements.empty() ) { // get all elements in the mesh SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator(); while ( eIt->more() ) - elems.insert( elems.end(), eIt->next()); + theElements.insert( theElements.end(), eIt->next()); } - else - elems = theElements; vector< TGroupOfElems > arrayOfGroups; TGroupOfElems groupOfElems; TMapOfNodeSet mapOfNodeSet; - TElemsSet::iterator elemIt = elems.begin(); - for ( int i = 0, j=0; elemIt != elems.end(); ++elemIt, ++j ) { + TIDSortedElemSet::iterator elemIt = theElements.begin(); + for ( int i = 0, j=0; elemIt != theElements.end(); ++elemIt, ++j ) { const SMDS_MeshElement* curElem = *elemIt; SortableElement SE(curElem); int ind = -1; @@ -7844,8 +8304,8 @@ void SMESH_MeshEditor::MergeElements(TListOfListOfElementsID & theGroupsOfElemen void SMESH_MeshEditor::MergeEqualElements() { - set aMeshElements; /* empty input - - to merge equal elements in the whole mesh */ + TIDSortedElemSet aMeshElements; /* empty input == + to merge equal elements in the whole mesh */ TListOfListOfElementsID aGroupsOfElementsID; FindEqualElements(aMeshElements, aGroupsOfElementsID); MergeElements(aGroupsOfElementsID); @@ -9225,7 +9685,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d) 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(); + aHelper.FixQuadraticElements(myError); } } @@ -9233,7 +9693,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d) /*! * \brief Makes given elements quadratic * \param theForce3d - if true, the medium nodes will be placed in the middle of link - * \param theElements - elements to make quadratic + * \param theElements - elements to make quadratic */ //================================================================================ @@ -9244,7 +9704,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, // we believe that all theElements are of the same type const SMDSAbs_ElementType elemType = (*theElements.begin())->GetType(); - + // get all nodes shared by theElements TIDSortedNodeSet allNodes; TIDSortedElemSet::iterator eIt = theElements.begin(); @@ -9365,7 +9825,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, if ( !theForce3d && !getenv("NO_FixQuadraticElements")) { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion helper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh - helper.FixQuadraticElements(); + helper.FixQuadraticElements( myError ); } } @@ -10229,7 +10689,7 @@ SMESH_MeshEditor::FindMatchingNodes(set& theSide1, \param theElems - the list of elements (edges or faces) to be replicated The nodes for duplication could be found from these elements \param theNodesNot - list of nodes to NOT replicate - \param theAffectedElems - the list of elements (cells and edges) to which the + \param theAffectedElems - the list of elements (cells and edges) to which the replicated nodes should be associated to. \return TRUE if operation has been completed successfully, FALSE otherwise */ @@ -10292,8 +10752,8 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh* theMeshDS, std::vector newNodes( anElem->NbNodes() ); SMDS_ElemIteratorPtr anIter = anElem->nodesIterator(); int ind = 0; - while ( anIter->more() ) - { + while ( anIter->more() ) + { SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next(); SMDS_MeshNode* aNewNode = aCurrNode; @@ -10328,14 +10788,14 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh* theMeshDS, /*! \brief Creates a hole in a mesh by doubling the nodes of some particular elements \param theNodes - identifiers of nodes to be doubled - \param theModifiedElems - identifiers of elements to be updated by the new (doubled) - nodes. If list of element identifiers is empty then nodes are doubled but + \param theModifiedElems - identifiers of elements to be updated by the new (doubled) + nodes. If list of element identifiers is empty then nodes are doubled but they not assigned to elements \return TRUE if operation has been completed successfully, FALSE otherwise */ //================================================================================ -bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, +bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, const std::list< int >& theListOfModifiedElems ) { MESSAGE("DoubleNodes"); @@ -10376,7 +10836,7 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, std::map< SMDS_MeshElement*, vector > anElemToNodes; std::list< int >::const_iterator anElemIter; - for ( anElemIter = theListOfModifiedElems.begin(); + for ( anElemIter = theListOfModifiedElems.begin(); anElemIter != theListOfModifiedElems.end(); ++anElemIter ) { int aCurr = *anElemIter; @@ -10388,8 +10848,8 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, SMDS_ElemIteratorPtr anIter = anElem->nodesIterator(); int ind = 0; - while ( anIter->more() ) - { + while ( anIter->more() ) + { SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next(); if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() ) { @@ -10402,7 +10862,7 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, anElemToNodes[ anElem ] = aNodeArr; } - // Change nodes of elements + // Change nodes of elements std::map< SMDS_MeshElement*, vector >::iterator anElemToNodesIter = anElemToNodes.begin(); @@ -10484,6 +10944,70 @@ namespace { }; } +//================================================================================ +/*! + \brief Identify the elements that will be affected by node duplication (actual duplication is not performed. + This method is the first step of DoubleNodeElemGroupsInRegion. + \param theElems - list of groups of elements (edges or faces) to be replicated + \param theNodesNot - list of groups of nodes not to replicated + \param theShape - shape to detect affected elements (element which geometric center + located on or inside shape). + The replicated nodes should be associated to affected elements. + \return groups of affected elements + \sa DoubleNodeElemGroupsInRegion() + */ +//================================================================================ + +bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems, + const TIDSortedElemSet& theNodesNot, + const TopoDS_Shape& theShape, + TIDSortedElemSet& theAffectedElems) +{ + if ( theShape.IsNull() ) + return false; + + const double aTol = Precision::Confusion(); + auto_ptr< BRepClass3d_SolidClassifier> bsc3d; + auto_ptr<_FaceClassifier> aFaceClassifier; + if ( theShape.ShapeType() == TopAbs_SOLID ) + { + bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));; + bsc3d->PerformInfinitePoint(aTol); + } + else if (theShape.ShapeType() == TopAbs_FACE ) + { + aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape))); + } + + // iterates on indicated elements and get elements by back references from their nodes + TIDSortedElemSet::const_iterator elemItr = theElems.begin(); + for ( ; elemItr != theElems.end(); ++elemItr ) + { + SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr; + if (!anElem) + continue; + + SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator(); + while ( nodeItr->more() ) + { + const SMDS_MeshNode* aNode = cast2Node(nodeItr->next()); + if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() ) + continue; + SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator(); + while ( backElemItr->more() ) + { + const SMDS_MeshElement* curElem = backElemItr->next(); + if ( curElem && theElems.find(curElem) == theElems.end() && + ( bsc3d.get() ? + isInside( curElem, *bsc3d, aTol ) : + isInside( curElem, *aFaceClassifier, aTol ))) + theAffectedElems.insert( curElem ); + } + } + } + return true; +} + //================================================================================ /*! \brief Creates a hole in a mesh by doubling the nodes of some particular elements @@ -10614,7 +11138,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector volume to modify) // with all the faces shared by 2 domains (group of elements) // and corresponding volume of this domain, for each shared face. - // a volume has a face shared by 2 domains if it has a neighbor which is not in is domain. + // a volume has a face shared by 2 domains if it has a neighbor which is not in his domain. //MESSAGE("Domain " << idom); const TIDSortedElemSet& domain = theElems[idom]; @@ -10677,8 +11201,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector cells; - cells.clear(); vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId); for (int i=0; i oldNodes; oldNodes.clear(); grid->GetNodeIds(oldNodes, face.cellId, face.cellType); - bool isMultipleDetected = false; std::set::iterator itn = oldNodes.begin(); for (; itn != oldNodes.end(); ++itn) { int oldId = *itn; - //MESSAGE(" node " << oldId); + //MESSAGE("-+-+-a node " << oldId); if (!nodeDomains.count(oldId)) nodeDomains[oldId] = emptyMap; // create an empty entry for node if (nodeDomains[oldId].empty()) - nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain + { + nodeDomains[oldId][idomain] = oldId; // keep the old node in the first domain + //MESSAGE("-+-+-b oldNode " << oldId << " domain " << idomain); + } std::map::iterator itdom = domvol.begin(); for (; itdom != domvol.end(); ++itdom) { @@ -10751,7 +11275,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector orderedDoms; //MESSAGE("multiple node " << oldId); - isMultipleDetected =true; if (mutipleNodes.count(oldId)) orderedDoms = mutipleNodes[oldId]; else @@ -10771,16 +11294,35 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorAddNode(coords[0], coords[1], coords[2]); int newId = newNode->getVtkId(); nodeDomains[oldId][idom] = newId; // cloned node for other domains - //MESSAGE(" newNode " << newId << " oldNode " << oldId << " size=" <= 3) - { - //MESSAGE("confirm multiple node " << oldId); - isMultipleDetected =true; + //MESSAGE("-+-+-c oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" < domvol = itface->second; + if (!domvol.count(idomain)) + continue; + DownIdType face = itface->first; + //MESSAGE(" --- face " << face.cellId); + std::set oldNodes; + oldNodes.clear(); + grid->GetNodeIds(oldNodes, face.cellId, face.cellType); + int nbMultipleNodes = 0; + std::set::iterator itn = oldNodes.begin(); + for (; itn != oldNodes.end(); ++itn) + { + int oldId = *itn; + if (mutipleNodes.count(oldId)) + nbMultipleNodes++; + } + if (nbMultipleNodes > 1) // check if an edge of the face is shared between 3 or more domains { //MESSAGE("multiple Nodes detected on a shared face"); int downId = itface->first.cellId; @@ -10794,7 +11336,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorgetDownArray(cellType)->getNumberOfDownCells(downId); @@ -10808,9 +11350,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector vn0 = mutipleNodes[nodes[0]]; vector vn1 = mutipleNodes[nodes[nbNodes - 1]]; - sort( vn0.begin(), vn0.end() ); - sort( vn1.begin(), vn1.end() ); - if (vn0 == vn1) + vector doms; + for (int i0 = 0; i0 < vn0.size(); i0++) + for (int i1 = 0; i1 < vn1.size(); i1++) + if (vn0[i0] == vn1[i1]) + doms.push_back(vn0[i0]); + if (doms.size() >2) { //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]); double *coords = grid->GetPoint(nodes[0]); @@ -10822,9 +11367,9 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector domvol; // domain --> a volume with the edge map angleDom; // oriented angles between planes defined by edge and volume centers int nbvol = grid->GetParentVolumes(vtkVolIds, downEdgeIds[ie], edgeType[ie]); - for (int id=0; id < vn0.size(); id++) + for (int id=0; id < doms.size(); id++) { - int idom = vn0[id]; + int idom = doms[id]; for (int ivol=0; ivolfromVtkToSmds(vtkVolIds[ivol]); @@ -10887,6 +11432,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectormyMesh->AddGroup(SMDSAbs_Face, joints2DName.c_str(), idg); + SMESHDS_Group *joints2DGrp = dynamic_cast(mapOfJunctionGroups[joints2DName]->GetGroupDS()); + string joints3DName = "joints3D"; + mapOfJunctionGroups[joints3DName] = this->myMesh->AddGroup(SMDSAbs_Volume, joints3DName.c_str(), idg); + SMESHDS_Group *joints3DGrp = dynamic_cast(mapOfJunctionGroups[joints3DName]->GetGroupDS()); + itface = faceDomains.begin(); for (; itface != faceDomains.end(); ++itface) { @@ -10910,13 +11463,16 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectormyMesh->AddGroup(vol->GetType(), namegrp.c_str(), idg); SMESHDS_Group *sgrp = dynamic_cast(mapOfJunctionGroups[namegrp]->GetGroupDS()); if (sgrp) sgrp->Add(vol->GetID()); + if (vol->GetType() == SMDSAbs_Volume) + joints3DGrp->Add(vol->GetID()); + else if (vol->GetType() == SMDSAbs_Face) + joints2DGrp->Add(vol->GetID()); } } @@ -10970,11 +11526,8 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorGetMeshDS()->AddVolumeFromVtkIds(orderedNodes); - stringstream grpname; - grpname << "mj_"; - grpname << 0 << "_" << 0; int idg; - string namegrp = grpname.str(); + string namegrp = "jointsMultiples"; if (!mapOfJunctionGroups.count(namegrp)) mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg); SMESHDS_Group *sgrp = dynamic_cast(mapOfJunctionGroups[namegrp]->GetGroupDS()); @@ -10983,7 +11536,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector& nodesCoords, + std::vector >& listOfListOfNodes) +{ + MESSAGE("--------------------------------"); + MESSAGE("SMESH_MeshEditor::CreateHoleSkin"); + MESSAGE("--------------------------------"); + + // --- zone of volumes to remove is given : + // 1 either by a geom shape (one or more vertices) and a radius, + // 2 either by a group of nodes (representative of the shape)to use with the radius, + // 3 either by a group of nodes where all the elements build on one of this nodes are to remove, + // In the case 2, the group of nodes is an external group of nodes from another mesh, + // In the case 3, the group of nodes is an internal group of the mesh (obtained for instance by a filter), + // defined by it's name. + + SMESHDS_GroupBase* groupDS = 0; + SMESH_Mesh::GroupIteratorPtr groupIt = this->myMesh->GetGroups(); + while ( groupIt->more() ) + { + groupDS = 0; + SMESH_Group * group = groupIt->next(); + if ( !group ) continue; + groupDS = group->GetGroupDS(); + if ( !groupDS || groupDS->IsEmpty() ) continue; + std::string grpName = group->GetName(); + //MESSAGE("grpName=" << grpName); + if (grpName == groupName) + break; + else + groupDS = 0; + } + + bool isNodeGroup = false; + bool isNodeCoords = false; + if (groupDS) + { + if (groupDS->GetType() != SMDSAbs_Node) + return; + isNodeGroup = true; // a group of nodes exists and it is in this mesh + } + + if (nodesCoords.size() > 0) + isNodeCoords = true; // a list o nodes given by their coordinates + //MESSAGE("---" << isNodeGroup << " " << isNodeCoords); + + // --- define groups to build + + int idg; // --- group of SMDS volumes + string grpvName = groupName; + grpvName += "_vol"; + SMESH_Group *grp = this->myMesh->AddGroup(SMDSAbs_Volume, grpvName.c_str(), idg); + if (!grp) + { + MESSAGE("group not created " << grpvName); + return; + } + SMESHDS_Group *sgrp = dynamic_cast(grp->GetGroupDS()); + + int idgs; // --- group of SMDS faces on the skin + string grpsName = groupName; + grpsName += "_skin"; + SMESH_Group *grps = this->myMesh->AddGroup(SMDSAbs_Face, grpsName.c_str(), idgs); + if (!grps) + { + MESSAGE("group not created " << grpsName); + return; + } + SMESHDS_Group *sgrps = dynamic_cast(grps->GetGroupDS()); + + int idgi; // --- group of SMDS faces internal (several shapes) + string grpiName = groupName; + grpiName += "_internalFaces"; + SMESH_Group *grpi = this->myMesh->AddGroup(SMDSAbs_Face, grpiName.c_str(), idgi); + if (!grpi) + { + MESSAGE("group not created " << grpiName); + return; + } + SMESHDS_Group *sgrpi = dynamic_cast(grpi->GetGroupDS()); + + int idgei; // --- group of SMDS faces internal (several shapes) + string grpeiName = groupName; + grpeiName += "_internalEdges"; + SMESH_Group *grpei = this->myMesh->AddGroup(SMDSAbs_Edge, grpeiName.c_str(), idgei); + if (!grpei) + { + MESSAGE("group not created " << grpeiName); + return; + } + SMESHDS_Group *sgrpei = dynamic_cast(grpei->GetGroupDS()); + + // --- build downward connectivity + + SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS(); + meshDS->BuildDownWardConnectivity(true); + SMDS_UnstructuredGrid* grid = meshDS->getGrid(); + + // --- set of volumes detected inside + + std::set setOfInsideVol; + std::set setOfVolToCheck; + + std::vector gpnts; + gpnts.clear(); + + if (isNodeGroup) // --- a group of nodes is provided : find all the volumes using one or more of this nodes + { + MESSAGE("group of nodes provided"); + SMDS_ElemIteratorPtr elemIt = groupDS->GetElements(); + while ( elemIt->more() ) + { + const SMDS_MeshElement* elem = elemIt->next(); + if (!elem) + continue; + const SMDS_MeshNode* node = dynamic_cast(elem); + if (!node) + continue; + SMDS_MeshElement* vol = 0; + SMDS_ElemIteratorPtr volItr = node->GetInverseElementIterator(SMDSAbs_Volume); + while (volItr->more()) + { + vol = (SMDS_MeshElement*)volItr->next(); + setOfInsideVol.insert(vol->getVtkId()); + sgrp->Add(vol->GetID()); + } + } + } + else if (isNodeCoords) + { + MESSAGE("list of nodes coordinates provided"); + int i = 0; + int k = 0; + while (i < nodesCoords.size()-2) + { + double x = nodesCoords[i++]; + double y = nodesCoords[i++]; + double z = nodesCoords[i++]; + gp_Pnt p = gp_Pnt(x, y ,z); + gpnts.push_back(p); + MESSAGE("TopoDS_Vertex " << k++ << " " << p.X() << " " << p.Y() << " " << p.Z()); + } + } + else // --- no group, no coordinates : use the vertices of the geom shape provided, and radius + { + MESSAGE("no group of nodes provided, using vertices from geom shape, and radius"); + TopTools_IndexedMapOfShape vertexMap; + TopExp::MapShapes( theShape, TopAbs_VERTEX, vertexMap ); + gp_Pnt p = gp_Pnt(0,0,0); + if (vertexMap.Extent() < 1) + return; + + for ( int i = 1; i <= vertexMap.Extent(); ++i ) + { + const TopoDS_Vertex& vertex = TopoDS::Vertex( vertexMap( i )); + p = BRep_Tool::Pnt(vertex); + gpnts.push_back(p); + MESSAGE("TopoDS_Vertex " << i << " " << p.X() << " " << p.Y() << " " << p.Z()); + } + } + + if (gpnts.size() > 0) + { + int nodeId = 0; + const SMDS_MeshNode* startNode = theNodeSearcher->FindClosestTo(gpnts[0]); + if (startNode) + nodeId = startNode->GetID(); + MESSAGE("nodeId " << nodeId); + + double radius2 = radius*radius; + MESSAGE("radius2 " << radius2); + + // --- volumes on start node + + setOfVolToCheck.clear(); + SMDS_MeshElement* startVol = 0; + SMDS_ElemIteratorPtr volItr = startNode->GetInverseElementIterator(SMDSAbs_Volume); + while (volItr->more()) + { + startVol = (SMDS_MeshElement*)volItr->next(); + setOfVolToCheck.insert(startVol->getVtkId()); + } + if (setOfVolToCheck.empty()) + { + MESSAGE("No volumes found"); + return; + } + + // --- starting with central volumes then their neighbors, check if they are inside + // or outside the domain, until no more new neighbor volume is inside. + // Fill the group of inside volumes + + std::map mapOfNodeDistance2; + mapOfNodeDistance2.clear(); + std::set setOfOutsideVol; + while (!setOfVolToCheck.empty()) + { + std::set::iterator it = setOfVolToCheck.begin(); + int vtkId = *it; + MESSAGE("volume to check, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + bool volInside = false; + vtkIdType npts = 0; + vtkIdType* pts = 0; + grid->GetCellPoints(vtkId, npts, pts); + for (int i=0; iGetPoint(pts[i]); + gp_Pnt aPoint = gp_Pnt(coords[0], coords[1], coords[2]); + distance2 = 1.E40; + for (int j=0; jAdd(meshDS->fromVtkToSmds(vtkId)); + break; + } + } + if (volInside) + { + setOfInsideVol.insert(vtkId); + MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + int neighborsVtkIds[NBMAXNEIGHBORS]; + int downIds[NBMAXNEIGHBORS]; + unsigned char downTypes[NBMAXNEIGHBORS]; + int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId); + for (int n = 0; n < nbNeighbors; n++) + if (!setOfInsideVol.count(neighborsVtkIds[n]) ||setOfOutsideVol.count(neighborsVtkIds[n])) + setOfVolToCheck.insert(neighborsVtkIds[n]); + } + else + { + setOfOutsideVol.insert(vtkId); + MESSAGE(" volume outside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + } + setOfVolToCheck.erase(vtkId); + } + } + + // --- for outside hexahedrons, check if they have more than one neighbor volume inside + // If yes, add the volume to the inside set + + bool addedInside = true; + std::set setOfVolToReCheck; + while (addedInside) + { + MESSAGE(" --------------------------- re check"); + addedInside = false; + std::set::iterator itv = setOfInsideVol.begin(); + for (; itv != setOfInsideVol.end(); ++itv) + { + int vtkId = *itv; + int neighborsVtkIds[NBMAXNEIGHBORS]; + int downIds[NBMAXNEIGHBORS]; + unsigned char downTypes[NBMAXNEIGHBORS]; + int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId); + for (int n = 0; n < nbNeighbors; n++) + if (!setOfInsideVol.count(neighborsVtkIds[n])) + setOfVolToReCheck.insert(neighborsVtkIds[n]); + } + setOfVolToCheck = setOfVolToReCheck; + setOfVolToReCheck.clear(); + while (!setOfVolToCheck.empty()) + { + std::set::iterator it = setOfVolToCheck.begin(); + int vtkId = *it; + if (grid->GetCellType(vtkId) == VTK_HEXAHEDRON) + { + MESSAGE("volume to recheck, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + int countInside = 0; + int neighborsVtkIds[NBMAXNEIGHBORS]; + int downIds[NBMAXNEIGHBORS]; + unsigned char downTypes[NBMAXNEIGHBORS]; + int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId); + for (int n = 0; n < nbNeighbors; n++) + if (setOfInsideVol.count(neighborsVtkIds[n])) + countInside++; + MESSAGE("countInside " << countInside); + if (countInside > 1) + { + MESSAGE(" volume inside, vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + setOfInsideVol.insert(vtkId); + sgrp->Add(meshDS->fromVtkToSmds(vtkId)); + addedInside = true; + } + else + setOfVolToReCheck.insert(vtkId); + } + setOfVolToCheck.erase(vtkId); + } + } + + // --- map of Downward faces at the boundary, inside the global volume + // map of Downward faces on the skin of the global volume (equivalent to SMDS faces on the skin) + // fill group of SMDS faces inside the volume (when several volume shapes) + // fill group of SMDS faces on the skin of the global volume (if skin) + + std::map boundaryFaces; // boundary faces inside the volume --> corresponding cell + std::map skinFaces; // faces on the skin of the global volume --> corresponding cell + std::set::iterator it = setOfInsideVol.begin(); + for (; it != setOfInsideVol.end(); ++it) + { + int vtkId = *it; + //MESSAGE(" vtkId " << vtkId << " smdsId " << meshDS->fromVtkToSmds(vtkId)); + int neighborsVtkIds[NBMAXNEIGHBORS]; + int downIds[NBMAXNEIGHBORS]; + unsigned char downTypes[NBMAXNEIGHBORS]; + int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId, true); + for (int n = 0; n < nbNeighbors; n++) + { + int neighborDim = SMDS_Downward::getCellDimension(grid->GetCellType(neighborsVtkIds[n])); + if (neighborDim == 3) + { + if (! setOfInsideVol.count(neighborsVtkIds[n])) // neighbor volume is not inside : face is boundary + { + DownIdType face(downIds[n], downTypes[n]); + boundaryFaces[face] = vtkId; + } + // if the face between to volumes is in the mesh, get it (internal face between shapes) + int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]); + if (vtkFaceId >= 0) + { + sgrpi->Add(meshDS->fromVtkToSmds(vtkFaceId)); + // find also the smds edges on this face + int nbEdges = grid->getDownArray(downTypes[n])->getNumberOfDownCells(downIds[n]); + const int* dEdges = grid->getDownArray(downTypes[n])->getDownCells(downIds[n]); + const unsigned char* dTypes = grid->getDownArray(downTypes[n])->getDownTypes(downIds[n]); + for (int i = 0; i < nbEdges; i++) + { + int vtkEdgeId = grid->getDownArray(dTypes[i])->getVtkCellId(dEdges[i]); + if (vtkEdgeId >= 0) + sgrpei->Add(meshDS->fromVtkToSmds(vtkEdgeId)); + } + } + } + else if (neighborDim == 2) // skin of the volume + { + DownIdType face(downIds[n], downTypes[n]); + skinFaces[face] = vtkId; + int vtkFaceId = grid->getDownArray(downTypes[n])->getVtkCellId(downIds[n]); + if (vtkFaceId >= 0) + sgrps->Add(meshDS->fromVtkToSmds(vtkFaceId)); + } + } + } + + // --- identify the edges constituting the wire of each subshape on the skin + // define polylines with the nodes of edges, equivalent to wires + // project polylines on subshapes, and partition, to get geom faces + + std::map > shapeIdToVtkIdSet; // shapeId --> set of vtkId on skin + std::set emptySet; + emptySet.clear(); + std::set shapeIds; + + SMDS_ElemIteratorPtr itelem = sgrps->GetElements(); + while (itelem->more()) + { + const SMDS_MeshElement *elem = itelem->next(); + int shapeId = elem->getshapeId(); + int vtkId = elem->getVtkId(); + if (!shapeIdToVtkIdSet.count(shapeId)) + { + shapeIdToVtkIdSet[shapeId] = emptySet; + shapeIds.insert(shapeId); + } + shapeIdToVtkIdSet[shapeId].insert(vtkId); + } + + std::map > shapeIdToEdges; // shapeId --> set of downward edges + std::set emptyEdges; + emptyEdges.clear(); + + std::map >::iterator itShape = shapeIdToVtkIdSet.begin(); + for (; itShape != shapeIdToVtkIdSet.end(); ++itShape) + { + int shapeId = itShape->first; + MESSAGE(" --- Shape ID --- "<< shapeId); + shapeIdToEdges[shapeId] = emptyEdges; + + std::vector nodesEdges; + + std::set::iterator its = itShape->second.begin(); + for (; its != itShape->second.end(); ++its) + { + int vtkId = *its; + MESSAGE(" " << vtkId); + int neighborsVtkIds[NBMAXNEIGHBORS]; + int downIds[NBMAXNEIGHBORS]; + unsigned char downTypes[NBMAXNEIGHBORS]; + int nbNeighbors = grid->GetNeighbors(neighborsVtkIds, downIds, downTypes, vtkId); + for (int n = 0; n < nbNeighbors; n++) + { + if (neighborsVtkIds[n]<0) // only smds faces are considered as neighbors here + continue; + int smdsId = meshDS->fromVtkToSmds(neighborsVtkIds[n]); + const SMDS_MeshElement* elem = meshDS->FindElement(smdsId); + if ( shapeIds.count(elem->getshapeId()) && !sgrps->Contains(elem)) // edge : neighbor in the set of shape, not in the group + { + DownIdType edge(downIds[n], downTypes[n]); + if (!shapeIdToEdges[shapeId].count(edge)) + { + shapeIdToEdges[shapeId].insert(edge); + int vtkNodeId[3]; + int nbNodes = grid->getDownArray(downTypes[n])->getNodes(downIds[n],vtkNodeId); + nodesEdges.push_back(vtkNodeId[0]); + nodesEdges.push_back(vtkNodeId[nbNodes-1]); + MESSAGE(" --- nodes " << vtkNodeId[0]+1 << " " << vtkNodeId[nbNodes-1]+1); + } + } + } + } + + std::list order; + order.clear(); + if (nodesEdges.size() > 0) + { + order.push_back(nodesEdges[0]); MESSAGE(" --- back " << order.back()+1); // SMDS id = VTK id + 1; + nodesEdges[0] = -1; + order.push_back(nodesEdges[1]); MESSAGE(" --- back " << order.back()+1); + nodesEdges[1] = -1; // do not reuse this edge + bool found = true; + while (found) + { + int nodeTofind = order.back(); // try first to push back + int i = 0; + for (i = 0; i use the previous one + if (nodesEdges[i-1] < 0) + found = false; + else + { + order.push_back(nodesEdges[i-1]); MESSAGE(" --- back " << order.back()+1); + nodesEdges[i-1] = -1; + } + else // even ==> use the next one + if (nodesEdges[i+1] < 0) + found = false; + else + { + order.push_back(nodesEdges[i+1]); MESSAGE(" --- back " << order.back()+1); + nodesEdges[i+1] = -1; + } + } + if (found) + continue; + // try to push front + found = true; + nodeTofind = order.front(); // try to push front + for (i = 0; i use the previous one + if (nodesEdges[i-1] < 0) + found = false; + else + { + order.push_front(nodesEdges[i-1]); MESSAGE(" --- front " << order.front()+1); + nodesEdges[i-1] = -1; + } + else // even ==> use the next one + if (nodesEdges[i+1] < 0) + found = false; + else + { + order.push_front(nodesEdges[i+1]); MESSAGE(" --- front " << order.front()+1); + nodesEdges[i+1] = -1; + } + } + } + + + std::vector nodes; + nodes.push_back(shapeId); + std::list::iterator itl = order.begin(); + for (; itl != order.end(); itl++) + { + nodes.push_back((*itl) + 1); // SMDS id = VTK id + 1; + MESSAGE(" ordered node " << nodes[nodes.size()-1]); + } + listOfListOfNodes.push_back(nodes); + } + + // partition geom faces with blocFissure + // mesh blocFissure and geom faces of the skin (external wires given, triangle algo to choose) + // mesh volume around blocFissure (skin triangles and quadrangle given, tetra algo to choose) + + return; +} + + //================================================================================ /*! * \brief Generates skin mesh (containing 2D cells) from 3D mesh @@ -11483,7 +12566,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, missType, /*noMedium=*/false)) continue; - SMDS_MeshElement* elem = + SMDS_MeshElement* elem = tgtEditor.AddElement(nodes, missType, !iQuad && nodes.size()/(iQuad+1)>4); ++nbAddedBnd; @@ -11533,7 +12616,7 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, { presentEditor->myLastCreatedElems.Append(presentBndElems[i]); } - + } // loop on given elements // ---------------------------------------------