X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=49f9c1b3a8776abf98730188553103d7e1ca84fb;hb=8c9a971309bd09be8c62a445303930dcff75ee3b;hp=b5c2af386a451fd74ead178ec2d79dbe01ba467f;hpb=45ed5d9a931bf24187e627f41e04baae23e7ffc4;p=modules%2Fsmesh.git diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index b5c2af386..49f9c1b3a 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -94,6 +94,8 @@ #include #include #include +#include +#include #define cast2Node(elem) static_cast( elem ) @@ -132,114 +134,125 @@ SMESH_MeshEditor::AddElement(const vector & node, int nbnode = node.size(); SMESHDS_Mesh* mesh = GetMeshDS(); switch ( type ) { - case SMDSAbs_0DElement: - if ( nbnode == 1 ) { - if ( ID >= 0 ) e = mesh->Add0DElementWithID(node[0], ID); - else e = mesh->Add0DElement (node[0] ); - } - break; - case SMDSAbs_Edge: - if ( nbnode == 2 ) { - if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], ID); - else e = mesh->AddEdge (node[0], node[1] ); - } - else if ( nbnode == 3 ) { - if ( ID >= 0 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID); - else e = mesh->AddEdge (node[0], node[1], node[2] ); - } - break; case SMDSAbs_Face: if ( !isPoly ) { if (nbnode == 3) { - if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID); - else e = mesh->AddFace (node[0], node[1], node[2] ); + if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], ID); + else e = mesh->AddFace (node[0], node[1], node[2] ); } else if (nbnode == 4) { - if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID); - else e = mesh->AddFace (node[0], node[1], node[2], node[3] ); + if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], ID); + else e = mesh->AddFace (node[0], node[1], node[2], node[3] ); } else if (nbnode == 6) { - if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], - node[4], node[5], ID); - else e = mesh->AddFace (node[0], node[1], node[2], node[3], - node[4], node[5] ); + if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], + node[4], node[5], ID); + else e = mesh->AddFace (node[0], node[1], node[2], node[3], + node[4], node[5] ); } else if (nbnode == 8) { - if ( ID >= 0 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], ID); - else e = mesh->AddFace (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7] ); + if ( ID >= 1 ) e = mesh->AddFaceWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], ID); + else e = mesh->AddFace (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7] ); } } else { - if ( ID >= 0 ) e = mesh->AddPolygonalFaceWithID(node, ID); - else e = mesh->AddPolygonalFace (node ); + if ( ID >= 1 ) e = mesh->AddPolygonalFaceWithID(node, ID); + else e = mesh->AddPolygonalFace (node ); } break; + case SMDSAbs_Volume: if ( !isPoly ) { if (nbnode == 4) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3] ); } else if (nbnode == 5) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4] ); } else if (nbnode == 6) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5] ); } else if (nbnode == 8) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7] ); } else if (nbnode == 10) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9] ); } else if (nbnode == 13) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12],ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12],ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12] ); } else if (nbnode == 15) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12],node[13],node[14],ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12],node[13],node[14] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12],node[13],node[14],ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12],node[13],node[14] ); } else if (nbnode == 20) { - if ( ID >= 0 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12],node[13],node[14],node[15], - node[16],node[17],node[18],node[19],ID); - else e = mesh->AddVolume (node[0], node[1], node[2], node[3], - node[4], node[5], node[6], node[7], - node[8], node[9], node[10],node[11], - node[12],node[13],node[14],node[15], - node[16],node[17],node[18],node[19] ); + if ( ID >= 1 ) e = mesh->AddVolumeWithID(node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12],node[13],node[14],node[15], + node[16],node[17],node[18],node[19],ID); + else e = mesh->AddVolume (node[0], node[1], node[2], node[3], + node[4], node[5], node[6], node[7], + node[8], node[9], node[10],node[11], + node[12],node[13],node[14],node[15], + node[16],node[17],node[18],node[19] ); } } + break; + + case SMDSAbs_Edge: + if ( nbnode == 2 ) { + if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], ID); + else e = mesh->AddEdge (node[0], node[1] ); + } + else if ( nbnode == 3 ) { + if ( ID >= 1 ) e = mesh->AddEdgeWithID(node[0], node[1], node[2], ID); + else e = mesh->AddEdge (node[0], node[1], node[2] ); + } + break; + + case SMDSAbs_0DElement: + if ( nbnode == 1 ) { + if ( ID >= 1 ) e = mesh->Add0DElementWithID(node[0], ID); + else e = mesh->Add0DElement (node[0] ); + } + break; + + case SMDSAbs_Node: + if ( ID >= 1 ) e = mesh->AddNodeWithID(node[0]->X(), node[0]->Y(), node[0]->Z(), ID); + else e = mesh->AddNode (node[0]->X(), node[0]->Y(), node[0]->Z()); + break; + + default:; } if ( e ) myLastCreatedElems.Append( e ); return e; @@ -352,45 +365,55 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem) if ( aMesh->ShapeToMesh().IsNull() ) return 0; - if ( theElem->GetType() == SMDSAbs_Node ) - { - int aShapeID = theElem->getshapeId(); - if (aShapeID <= 0) - return 0; - else - return aShapeID; - } - - TopoDS_Shape aShape; // the shape a node is on - SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); - while ( nodeIt->more() ) { - const SMDS_MeshNode* node = static_cast( nodeIt->next() ); - int aShapeID = node->getshapeId(); - if (aShapeID > 0) { - SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ); - if ( sm ) { - if ( sm->Contains( theElem )) - return aShapeID; - if ( aShape.IsNull() ) - aShape = aMesh->IndexToShape( aShapeID ); - } - else { - //MESSAGE ( "::FindShape() No SubShape for aShapeID " << aShapeID ); + int aShapeID = theElem->getshapeId(); + if ( aShapeID < 1 ) + return 0; + + if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID )) + if ( sm->Contains( theElem )) + return aShapeID; + + if ( theElem->GetType() == SMDSAbs_Node ) { + MESSAGE( ":( Error: invalid myShapeId of node " << theElem->GetID() ); + } + else { + MESSAGE( ":( Error: invalid myShapeId of element " << theElem->GetID() ); + } + + TopoDS_Shape aShape; // the shape a node of theElem is on + if ( theElem->GetType() != SMDSAbs_Node ) + { + SMDS_ElemIteratorPtr nodeIt = theElem->nodesIterator(); + while ( nodeIt->more() ) { + const SMDS_MeshNode* node = static_cast( nodeIt->next() ); + if ((aShapeID = node->getshapeId()) > 0) { + if ( SMESHDS_SubMesh * sm = aMesh->MeshElements( aShapeID ) ) { + if ( sm->Contains( theElem )) + return aShapeID; + if ( aShape.IsNull() ) + aShape = aMesh->IndexToShape( aShapeID ); + } } } } // None of nodes is on a proper shape, // find the shape among ancestors of aShape on which a node is - if ( aShape.IsNull() ) { - //MESSAGE ("::FindShape() - NONE node is on shape") - return 0; + if ( !aShape.IsNull() ) { + TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape )); + for ( ; ancIt.More(); ancIt.Next() ) { + SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() ); + if ( sm && sm->Contains( theElem )) + return aMesh->ShapeToIndex( ancIt.Value() ); + } } - TopTools_ListIteratorOfListOfShape ancIt( GetMesh()->GetAncestors( aShape )); - for ( ; ancIt.More(); ancIt.Next() ) { - SMESHDS_SubMesh * sm = aMesh->MeshElements( ancIt.Value() ); - if ( sm && sm->Contains( theElem )) - return aMesh->ShapeToIndex( ancIt.Value() ); + else + { + const map& id2sm = GetMeshDS()->SubMeshes(); + map::const_iterator id_sm = id2sm.begin(); + for ( ; id_sm != id2sm.end(); ++id_sm ) + if ( id_sm->second->Contains( theElem )) + return id_sm->first; } //MESSAGE ("::FindShape() - SHAPE NOT FOUND") @@ -1552,14 +1575,14 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, const int theMethodFlags) { // std-like iterator on coordinates of nodes of mesh element - typedef SMDS_StdIterator< TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator; + typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator; NXyzIterator xyzEnd; SMDS_VolumeTool volTool; SMESH_MesherHelper helper( *GetMesh()); - SMESHDS_SubMesh* subMesh = GetMeshDS()->MeshElements(1); - SMESHDS_SubMesh* fSubMesh = subMesh; + SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1); + SMESHDS_SubMesh* fSubMesh = 0;//subMesh; SMESH_SequenceOfElemPtr newNodes, newElems; @@ -1717,10 +1740,10 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } for ( int i = 0; i < triangles.size(); ++i ) { - if ( !triangles.back() ) continue; + if ( !triangles[i] ) continue; if ( fSubMesh ) - fSubMesh->AddElement( triangles.back()); - newElems.Append( triangles.back() ); + fSubMesh->AddElement( triangles[i]); + newElems.Append( triangles[i] ); } ReplaceElemInGroups( face, triangles, GetMeshDS() ); GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false ); @@ -2670,6 +2693,9 @@ void SMESH_MeshEditor::GetLinkedNodes( const SMDS_MeshNode* theNode, while ( elemIt->more() ) { const SMDS_MeshElement* elem = elemIt->next(); + if(elem->GetType() == SMDSAbs_0DElement) + continue; + SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); if ( elem->GetType() == SMDSAbs_Volume ) { @@ -3874,6 +3900,8 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, const SMDS_MeshElement* elem = itElem->first; vector& vecNewNodes = itElemNodes->second; + if(itElem->second.size()==0) continue; + if ( elem->GetType() == SMDSAbs_Edge ) { // create a ceiling edge if (!elem->IsQuadratic()) { @@ -3898,8 +3926,6 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, if ( elem->GetType() != SMDSAbs_Face ) continue; - if(itElem->second.size()==0) continue; - bool hasFreeLinks = false; TIDSortedElemSet avoidSet; @@ -5292,7 +5318,7 @@ void SMESH_MeshEditor::LinearAngleVariation(const int nbSteps, * \param theCopy - if true, create translated copies of theElems * \param theMakeGroups - if true and theCopy, create translated groups * \param theTargetMesh - mesh to copy translated elements into - * \retval SMESH_MeshEditor::PGroupIDs - list of ids of created groups + * \return SMESH_MeshEditor::PGroupIDs - list of ids of created groups */ //================================================================================ @@ -5310,22 +5336,37 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems, string groupPostfix; switch ( theTrsf.Form() ) { case gp_PntMirror: + MESSAGE("gp_PntMirror"); + needReverse = true; + groupPostfix = "mirrored"; + break; case gp_Ax1Mirror: + MESSAGE("gp_Ax1Mirror"); + groupPostfix = "mirrored"; + break; case gp_Ax2Mirror: + MESSAGE("gp_Ax2Mirror"); needReverse = true; groupPostfix = "mirrored"; break; case gp_Rotation: + MESSAGE("gp_Rotation"); groupPostfix = "rotated"; break; case gp_Translation: + MESSAGE("gp_Translation"); groupPostfix = "translated"; break; case gp_Scale: + MESSAGE("gp_Scale"); + groupPostfix = "scaled"; + break; case gp_CompoundTrsf: // different scale by axis + MESSAGE("gp_CompoundTrsf"); groupPostfix = "scaled"; break; default: + MESSAGE("default"); needReverse = false; groupPostfix = "transformed"; } @@ -6215,7 +6256,7 @@ struct SMESH_NodeSearcherImpl: public SMESH_NodeSearcher const SMDS_MeshNode* closestNode = 0; list::iterator nIt = nodes.begin(); for ( ; nIt != nodes.end(); ++nIt ) { - double sqDist = thePnt.SquareDistance( SMESH_MeshEditor::TNodeXYZ( *nIt ) ); + double sqDist = thePnt.SquareDistance( SMESH_TNodeXYZ( *nIt ) ); if ( minSqDist > sqDist ) { closestNode = *nIt; minSqDist = sqDist; @@ -6270,7 +6311,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() { public: - ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance = NodeRadius ); + 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); ~ElementBndBoxTree(); @@ -6297,13 +6338,13 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() */ //================================================================================ - ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, double tolerance) + ElementBndBoxTree::ElementBndBoxTree(const SMDS_Mesh& mesh, SMDSAbs_ElementType elemType, SMDS_ElemIteratorPtr theElemIt, double tolerance) :SMESH_Octree( new SMESH_Octree::Limit( MaxLevel, /*minSize=*/0. )) { int nbElems = mesh.GetMeshInfo().NbElements( elemType ); _elements.reserve( nbElems ); - SMDS_ElemIteratorPtr elemIt = mesh.elementsIterator( elemType ); + SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); while ( elemIt->more() ) _elements.push_back( new ElementBox( elemIt->next(),tolerance )); @@ -6435,7 +6476,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() _refCount = 1; SMDS_ElemIteratorPtr nIt = elem->nodesIterator(); while ( nIt->more() ) - Add( SMESH_MeshEditor::TNodeXYZ( cast2Node( nIt->next() ))); + Add( SMESH_TNodeXYZ( cast2Node( nIt->next() ))); Enlarge( tolerance ); } @@ -6451,6 +6492,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher { SMESHDS_Mesh* _mesh; + SMDS_ElemIteratorPtr _meshPartIt; ElementBndBoxTree* _ebbTree; SMESH_NodeSearcherImpl* _nodeSearcher; SMDSAbs_ElementType _elementType; @@ -6458,8 +6500,8 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher bool _outerFacesFound; set _outerFaces; // empty means "no internal faces at all" - SMESH_ElementSearcherImpl( SMESHDS_Mesh& mesh ) - : _mesh(&mesh),_ebbTree(0),_nodeSearcher(0), _tolerance(-1), _outerFacesFound(false) {} + 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; @@ -6541,7 +6583,7 @@ double SMESH_ElementSearcherImpl::getTolerance() SMDS_NodeIteratorPtr nodeIt = _mesh->nodesIterator(); elemSize = 1; if ( meshInfo.NbNodes() > 2 ) - elemSize = SMESH_MeshEditor::TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() ); + elemSize = SMESH_TNodeXYZ( nodeIt->next() ).Distance( nodeIt->next() ); } else { @@ -6549,7 +6591,8 @@ double SMESH_ElementSearcherImpl::getTolerance() _mesh->elementsIterator( SMDSAbs_ElementType( complexType )); const SMDS_MeshElement* elem = elemIt->next(); SMDS_ElemIteratorPtr nodeIt = elem->nodesIterator(); - SMESH_MeshEditor::TNodeXYZ n1( cast2Node( nodeIt->next() )); + SMESH_TNodeXYZ n1( cast2Node( nodeIt->next() )); + elemSize = 0; while ( nodeIt->more() ) { double dist = n1.Distance( cast2Node( nodeIt->next() )); @@ -6582,8 +6625,8 @@ bool SMESH_ElementSearcherImpl::getIntersParamOnLine(const gp_Lin& lin int nbNodes = face->IsQuadratic() ? face->NbNodes()/2 : face->NbNodes(); for ( int i = 0; i < nbNodes && nbInts < 2; ++i ) { - GC_MakeSegment edge( SMESH_MeshEditor::TNodeXYZ( face->GetNode( i )), - SMESH_MeshEditor::TNodeXYZ( face->GetNode( (i+1)%nbNodes) )); + 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) { @@ -6643,8 +6686,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF seamLinks.insert( link ); // link direction within the outerFace - gp_Vec n1n2( SMESH_MeshEditor::TNodeXYZ( link.node1()), - SMESH_MeshEditor::TNodeXYZ( link.node2())); + 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 ); @@ -6727,7 +6770,7 @@ FindElementsByPoint(const gp_Pnt& point, const SMDS_MeshNode* closeNode = _nodeSearcher->FindClosestTo( point ); if ( !closeNode ) return foundElements.size(); - if ( point.Distance( SMESH_MeshEditor::TNodeXYZ( closeNode )) > tolerance ) + if ( point.Distance( SMESH_TNodeXYZ( closeNode )) > tolerance ) return foundElements.size(); // to far from any node if ( type == SMDSAbs_Node ) @@ -6747,7 +6790,7 @@ FindElementsByPoint(const gp_Pnt& point, if ( !_ebbTree || _elementType != type ) { if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, tolerance ); + _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt, tolerance ); } TIDSortedElemSet suspectElems; _ebbTree->getElementsNearPoint( point, suspectElems ); @@ -6771,7 +6814,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) if ( !_ebbTree || _elementType != SMDSAbs_Face ) { if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = SMDSAbs_Face ); + _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 @@ -6799,7 +6842,7 @@ TopAbs_State SMESH_ElementSearcherImpl::GetPointState(const gp_Pnt& point) // get face plane gp_XYZ fNorm; if ( !SMESH_Algo::FaceNormal( *face, fNorm, /*normalized=*/false)) continue; - gp_Pln facePlane( SMESH_MeshEditor::TNodeXYZ( (*face)->GetNode(0)), fNorm ); + gp_Pln facePlane( SMESH_TNodeXYZ( (*face)->GetNode(0)), fNorm ); // perform intersection IntAna_IntConicQuad intersection( line, IntAna_Quadric( facePlane )); @@ -6995,7 +7038,7 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& if ( !_ebbTree || _elementType != type ) { if ( _ebbTree ) delete _ebbTree; - _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type ); + _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt ); } TIDSortedElemSet suspectFaces; // elements possibly intersecting the line _ebbTree->getElementsNearLine( line, suspectFaces ); @@ -7013,6 +7056,17 @@ SMESH_ElementSearcher* SMESH_MeshEditor::GetElementSearcher() return new SMESH_ElementSearcherImpl( *GetMeshDS() ); } +//======================================================================= +/*! + * \brief Return SMESH_ElementSearcher + */ +//======================================================================= + +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 @@ -7041,7 +7095,7 @@ bool SMESH_MeshEditor::isOut( const SMDS_MeshElement* element, const gp_Pnt& poi while ( nodeIt->more() ) { const SMDS_MeshNode* node = cast2Node( nodeIt->next() ); - xyz.push_back( TNodeXYZ(node) ); + xyz.push_back( SMESH_TNodeXYZ(node) ); nodeList.push_back(node); } @@ -7368,10 +7422,9 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) vector polygons_nodes; vector quantities; int nbNew = SimplifyFace(face_nodes, polygons_nodes, quantities); - if (nbNew > 0) { inode = 0; - for (int iface = 0; iface < nbNew - 1; iface++) { + for (int iface = 0; iface < nbNew; iface++) { int nbNodes = quantities[iface]; vector poly_nodes (nbNodes); for (int ii = 0; ii < nbNodes; ii++, inode++) { @@ -7838,7 +7891,11 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes) //MESSAGE("Change regular element or polygon " << elemId); SMDSAbs_ElementType etyp = elem->GetType(); uniqueNodes.resize(nbUniqueNodes); - SMDS_MeshElement* newElem = this->AddElement(uniqueNodes, etyp, false); + SMDS_MeshElement* newElem = 0; + if (elem->GetEntityType() == SMDSEntity_Polygon) + newElem = this->AddElement(uniqueNodes, etyp, true); + else + newElem = this->AddElement(uniqueNodes, etyp, false); if (newElem) { myLastCreatedElems.Append(newElem); @@ -9138,7 +9195,7 @@ void SMESH_MeshEditor::UpdateVolumes (const SMDS_MeshNode* theBetweenNode //======================================================================= /*! * \brief Convert elements contained in a submesh to quadratic - * \retval int - nb of checked elements + * \return int - nb of checked elements */ //======================================================================= @@ -9158,8 +9215,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, if( !elem || elem->IsQuadratic() ) continue; int id = elem->GetID(); - //MESSAGE("elem " << id); - id = 0; // get a free number for new elements int nbNodes = elem->NbNodes(); SMDSAbs_ElementType aType = elem->GetType(); @@ -9167,6 +9222,8 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, if ( elem->GetEntityType() == SMDSEntity_Polyhedra ) nbNodeInFaces = static_cast( elem )->GetQuantities(); + GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false); + const SMDS_MeshElement* NewElem = 0; switch( aType ) @@ -9220,8 +9277,6 @@ int SMESH_MeshEditor::convertElemToQuadratic(SMESHDS_SubMesh * theSm, ReplaceElemInGroups( elem, NewElem, GetMeshDS()); if( NewElem ) theSm->AddElement( NewElem ); - - GetMeshDS()->RemoveFreeElement(elem, theSm, /*fromGroups=*/false); } // if (!GetMeshDS()->isCompacted()) // GetMeshDS()->compactMesh(); @@ -9347,14 +9402,152 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d) aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh aHelper.FixQuadraticElements(); } - if (!GetMeshDS()->isCompacted()) - GetMeshDS()->compactMesh(); +} + +//================================================================================ +/*! + * \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 + */ +//================================================================================ + +void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, + TIDSortedElemSet& theElements) +{ + if ( theElements.empty() ) return; + + // we believe that all theElements are of the same type + SMDSAbs_ElementType elemType = (*theElements.begin())->GetType(); + + // get all nodes shared by theElements + TIDSortedNodeSet allNodes; + TIDSortedElemSet::iterator eIt = theElements.begin(); + for ( ; eIt != theElements.end(); ++eIt ) + allNodes.insert( (*eIt)->begin_nodes(), (*eIt)->end_nodes() ); + + // complete theElements with elements of lower dim whose all nodes are in allNodes + + TIDSortedElemSet quadAdjacentElems [ SMDSAbs_NbElementTypes ]; // quadratic adjacent elements + TIDSortedElemSet checkedAdjacentElems [ SMDSAbs_NbElementTypes ]; + TIDSortedNodeSet::iterator nIt = allNodes.begin(); + for ( ; nIt != allNodes.end(); ++nIt ) + { + const SMDS_MeshNode* n = *nIt; + SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator(); + while ( invIt->more() ) + { + const SMDS_MeshElement* e = invIt->next(); + if ( e->IsQuadratic() ) + { + quadAdjacentElems[ e->GetType() ].insert( e ); + continue; + } + if ( e->GetType() >= elemType ) + { + continue; // same type of more complex linear element + } + + if ( !checkedAdjacentElems[ e->GetType() ].insert( e ).second ) + continue; // e is already checked + + // check nodes + bool allIn = true; + SMDS_ElemIteratorPtr nodeIt = e->nodesIterator(); + while ( nodeIt->more() && allIn ) + allIn = allNodes.count( cast2Node( nodeIt->next() )); + if ( allIn ) + theElements.insert(e ); + } + } + + SMESH_MesherHelper helper(*myMesh); + helper.SetIsQuadratic( true ); + + // add links of quadratic adjacent elements to the helper + + if ( !quadAdjacentElems[SMDSAbs_Edge].empty() ) + for ( eIt = quadAdjacentElems[SMDSAbs_Edge].begin(); + eIt != quadAdjacentElems[SMDSAbs_Edge].end(); ++eIt ) + { + helper.AddTLinks( static_cast< const SMDS_MeshEdge*> (*eIt) ); + } + if ( !quadAdjacentElems[SMDSAbs_Face].empty() ) + for ( eIt = quadAdjacentElems[SMDSAbs_Face].begin(); + eIt != quadAdjacentElems[SMDSAbs_Face].end(); ++eIt ) + { + helper.AddTLinks( static_cast< const SMDS_MeshFace*> (*eIt) ); + } + if ( !quadAdjacentElems[SMDSAbs_Volume].empty() ) + for ( eIt = quadAdjacentElems[SMDSAbs_Volume].begin(); + eIt != quadAdjacentElems[SMDSAbs_Volume].end(); ++eIt ) + { + helper.AddTLinks( static_cast< const SMDS_MeshVolume*> (*eIt) ); + } + + // make quadratic elements instead of linear ones + + SMESHDS_Mesh* meshDS = GetMeshDS(); + SMESHDS_SubMesh* smDS = 0; + for ( eIt = theElements.begin(); eIt != theElements.end(); ++eIt ) + { + const SMDS_MeshElement* elem = *eIt; + if( elem->IsQuadratic() || elem->NbNodes() < 2 || elem->IsPoly() ) + continue; + + int id = elem->GetID(); + SMDSAbs_ElementType type = elem->GetType(); + vector nodes ( elem->begin_nodes(), elem->end_nodes()); + + if ( !smDS || !smDS->Contains( elem )) + smDS = meshDS->MeshElements( elem->getshapeId() ); + meshDS->RemoveFreeElement(elem, smDS, /*fromGroups=*/false); + + SMDS_MeshElement * newElem = 0; + switch( nodes.size() ) + { + case 4: // cases for most multiple element types go first (for optimization) + if ( type == SMDSAbs_Volume ) + newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d); + else + newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], nodes[3], id, theForce3d); + break; + case 8: + newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5], nodes[6], nodes[7], id, theForce3d); + break; + case 3: + newElem = helper.AddFace (nodes[0], nodes[1], nodes[2], id, theForce3d); + break; + case 2: + newElem = helper.AddEdge(nodes[0], nodes[1], id, theForce3d); + break; + case 5: + newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], id, theForce3d); + break; + case 6: + newElem = helper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], + nodes[4], nodes[5], id, theForce3d); + break; + default:; + } + ReplaceElemInGroups( elem, newElem, meshDS); + if( newElem && smDS ) + smDS->AddElement( newElem ); + } + + 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(); + } } //======================================================================= /*! * \brief Convert quadratic elements to linear ones and remove quadratic nodes - * \retval int - nb of checked elements + * \return int - nb of checked elements */ //======================================================================= @@ -9364,7 +9557,6 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm, { int nbElem = 0; SMESHDS_Mesh* meshDS = GetMeshDS(); - const bool notFromGroups = false; while( theItr->more() ) { @@ -9372,44 +9564,28 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm, nbElem++; if( elem && elem->IsQuadratic()) { - int id = elem->GetID(); - int nbNodes = elem->NbNodes(); - vector nodes, mediumNodes; - nodes.reserve( nbNodes ); - mediumNodes.reserve( nbNodes ); - - for(int i = 0; i < nbNodes; i++) - { - const SMDS_MeshNode* n = elem->GetNode(i); - - if( elem->IsMediumNode( n ) ) - mediumNodes.push_back( n ); - else - nodes.push_back( n ); - } - if( nodes.empty() ) continue; + int id = elem->GetID(); + int nbCornerNodes = elem->NbCornerNodes(); SMDSAbs_ElementType aType = elem->GetType(); - //remove old quadratic element - meshDS->RemoveFreeElement( elem, theSm, notFromGroups ); + vector nodes( elem->begin_nodes(), elem->end_nodes() ); - SMDS_MeshElement * NewElem = AddElement( nodes, aType, false, id ); - ReplaceElemInGroups(elem, NewElem, meshDS); - if( theSm && NewElem ) - theSm->AddElement( NewElem ); + //remove a quadratic element + if ( !theSm || !theSm->Contains( elem )) + theSm = meshDS->MeshElements( elem->getshapeId() ); + meshDS->RemoveFreeElement( elem, theSm, /*fromGroups=*/false ); // remove medium nodes - vector::iterator nIt = mediumNodes.begin(); - for ( ; nIt != mediumNodes.end(); ++nIt ) { - const SMDS_MeshNode* n = *nIt; - if ( n->NbInverseElements() == 0 ) { - if ( n->getshapeId() != theShapeID ) - meshDS->RemoveFreeNode( n, meshDS->MeshElements - ( n->getshapeId() )); - else - meshDS->RemoveFreeNode( n, theSm ); - } - } + for ( unsigned i = nbCornerNodes; i < nodes.size(); ++i ) + if ( nodes[i]->NbInverseElements() == 0 ) + meshDS->RemoveFreeNode( nodes[i], theSm ); + + // add a linear element + nodes.resize( nbCornerNodes ); + SMDS_MeshElement * newElem = AddElement( nodes, aType, false, id ); + ReplaceElemInGroups(elem, newElem, meshDS); + if( theSm && newElem ) + theSm->AddElement( newElem ); } } return nbElem; @@ -9419,7 +9595,8 @@ int SMESH_MeshEditor::removeQuadElem(SMESHDS_SubMesh * theSm, //function : ConvertFromQuadratic //purpose : //======================================================================= -bool SMESH_MeshEditor::ConvertFromQuadratic() + +bool SMESH_MeshEditor::ConvertFromQuadratic() { int nbCheckedElems = 0; if ( myMesh->HasShapeToMesh() ) @@ -9446,6 +9623,102 @@ bool SMESH_MeshEditor::ConvertFromQuadratic() return true; } +namespace +{ + //================================================================================ + /*! + * \brief Return true if all medium nodes of the element are in the node set + */ + //================================================================================ + + bool allMediumNodesIn(const SMDS_MeshElement* elem, TIDSortedNodeSet& nodeSet ) + { + for ( int i = elem->NbCornerNodes(); i < elem->NbNodes(); ++i ) + if ( !nodeSet.count( elem->GetNode(i) )) + return false; + return true; + } +} + +//================================================================================ +/*! + * \brief Makes given elements linear + */ +//================================================================================ + +void SMESH_MeshEditor::ConvertFromQuadratic(TIDSortedElemSet& theElements) +{ + if ( theElements.empty() ) return; + + // collect IDs of medium nodes of theElements; some of these nodes will be removed + set mediumNodeIDs; + TIDSortedElemSet::iterator eIt = theElements.begin(); + for ( ; eIt != theElements.end(); ++eIt ) + { + const SMDS_MeshElement* e = *eIt; + for ( int i = e->NbCornerNodes(); i < e->NbNodes(); ++i ) + mediumNodeIDs.insert( e->GetNode(i)->GetID() ); + } + + // replace given elements by linear ones + typedef SMDS_SetIterator TSetIterator; + SMDS_ElemIteratorPtr elemIt( new TSetIterator( theElements.begin(), theElements.end() )); + removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 ); + + // we need to convert remaining elements whose all medium nodes are in mediumNodeIDs + // except those elements sharing medium nodes of quadratic element whose medium nodes + // are not all in mediumNodeIDs + + // get remaining medium nodes + TIDSortedNodeSet mediumNodes; + set::iterator nIdsIt = mediumNodeIDs.begin(); + for ( ; nIdsIt != mediumNodeIDs.end(); ++nIdsIt ) + if ( const SMDS_MeshNode* n = GetMeshDS()->FindNode( *nIdsIt )) + mediumNodes.insert( mediumNodes.end(), n ); + + // find more quadratic elements to convert + TIDSortedElemSet moreElemsToConvert; + TIDSortedNodeSet::iterator nIt = mediumNodes.begin(); + for ( ; nIt != mediumNodes.end(); ++nIt ) + { + SMDS_ElemIteratorPtr invIt = (*nIt)->GetInverseElementIterator(); + while ( invIt->more() ) + { + const SMDS_MeshElement* e = invIt->next(); + if ( e->IsQuadratic() && allMediumNodesIn( e, mediumNodes )) + { + // find a more complex element including e and + // whose medium nodes are not in mediumNodes + bool complexFound = false; + for ( int type = e->GetType() + 1; type < SMDSAbs_0DElement; ++type ) + { + SMDS_ElemIteratorPtr invIt2 = + (*nIt)->GetInverseElementIterator( SMDSAbs_ElementType( type )); + while ( invIt2->more() ) + { + const SMDS_MeshElement* eComplex = invIt2->next(); + if ( eComplex->IsQuadratic() && !allMediumNodesIn( eComplex, mediumNodes)) + { + int nbCommonNodes = SMESH_Algo::GetCommonNodes( e, eComplex ).size(); + if ( nbCommonNodes == e->NbNodes()) + { + complexFound = true; + type = SMDSAbs_NbElementTypes; // to quit from the outer loop + break; + } + } + } + } + if ( !complexFound ) + moreElemsToConvert.insert( e ); + } + } + } + elemIt = SMDS_ElemIteratorPtr + (new TSetIterator( moreElemsToConvert.begin(), moreElemsToConvert.end() )); + removeQuadElem( /*theSm=*/0, elemIt, /*theShapeID=*/0 ); +} + //======================================================================= //function : SewSideElements //purpose : @@ -10075,7 +10348,7 @@ SMESH_MeshEditor::SewSideElements (TIDSortedElemSet& theSide1, * \param theSecondNode1 - a boundary node of set 1 linked with theFirstNode1 * \param theSecondNode2 - a node of set 2 corresponding to theSecondNode1 * \param nReplaceMap - output map of corresponding nodes - * \retval bool - is a success or not + * \return bool - is a success or not */ //================================================================================ @@ -10452,7 +10725,7 @@ namespace { gp_XYZ centerXYZ (0, 0, 0); SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator(); while (aNodeItr->more()) - centerXYZ += SMESH_MeshEditor::TNodeXYZ(cast2Node( aNodeItr->next())); + centerXYZ += SMESH_TNodeXYZ(cast2Node( aNodeItr->next())); gp_Pnt aPnt = centerXYZ / theElem->NbNodes(); theClassifier.Perform(aPnt, theTol); @@ -10557,12 +10830,35 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems, return DoubleNodes( theElems, theNodesNot, anAffected ); } +/*! + * \brief compute an oriented angle between two planes defined by four points. + * The vector (p0,p1) defines the intersection of the 2 planes (p0,p1,g1) and (p0,p1,g2) + * @param p0 base of the rotation axe + * @param p1 extremity of the rotation axe + * @param g1 belongs to the first plane + * @param g2 belongs to the second plane + */ +double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Pnt& g1, const gp_Pnt& g2) +{ +// MESSAGE(" p0: " << p0.X() << " " << p0.Y() << " " << p0.Z()); +// MESSAGE(" p1: " << p1.X() << " " << p1.Y() << " " << p1.Z()); +// MESSAGE(" g1: " << g1.X() << " " << g1.Y() << " " << g1.Z()); +// MESSAGE(" g2: " << g2.X() << " " << g2.Y() << " " << g2.Z()); + gp_Vec vref(p0, p1); + gp_Vec v1(p0, g1); + gp_Vec v2(p0, g2); + gp_Vec n1 = vref.Crossed(v1); + gp_Vec n2 = vref.Crossed(v2); + return n2.AngleWithRef(n1, vref); +} + /*! * \brief Double nodes on shared faces between groups of volumes and create flat elements on demand. * The list of groups must describe a partition of the mesh volumes. * The nodes of the internal faces at the boundaries of the groups are doubled. * In option, the internal faces are replaced by flat elements. * Triangles are transformed in prisms, and quadrangles in hexahedrons. + * The flat elements are stored in groups of volumes. * @param theElems - list of groups of volumes, where a group of volume is a set of * SMDS_MeshElements sorted by Id. * @param createJointElems - if TRUE, create the elements @@ -10571,23 +10867,29 @@ bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems, bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector& theElems, bool createJointElems) { - MESSAGE("------------------------------------------------------"); - MESSAGE("SMESH_MeshEditor::CreateJointElementsOnGroupBoundaries"); - MESSAGE("------------------------------------------------------"); + MESSAGE("----------------------------------------------"); + MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries"); + MESSAGE("----------------------------------------------"); SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS(); - meshDS->BuildDownWardConnectivity(false); + meshDS->BuildDownWardConnectivity(true); CHRONO(50); SMDS_UnstructuredGrid *grid = meshDS->getGrid(); // --- build the list of faces shared by 2 domains (group of elements), with their domain and volume indexes + // build the list of cells with only a node or an edge on the border, with their domain and volume indexes // build the list of nodes shared by 2 or more domains, with their domain indexes - std::map, DownIdCompare> faceDomains; // 2x(id domain --> id volume) - std::map > nodeDomains; //oldId -> (domainId -> newId) + std::map, DownIdCompare> faceDomains; // face --> (id domain --> id volume) + std::mapcelldom; // cell vtkId --> domain + std::map, DownIdCompare> cellDomains; // oldNode --> (id domain --> id cell) + std::map > nodeDomains; // oldId --> (domainId --> newId) faceDomains.clear(); + celldom.clear(); + cellDomains.clear(); nodeDomains.clear(); std::map emptyMap; + std::set emptySet; emptyMap.clear(); for (int idom = 0; idom < theElems.size(); idom++) @@ -10622,43 +10924,217 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector, DownIdCompare>::iterator itface; + + // --- explore the shared faces domain by domain, + // explore the nodes of the face and see if they belong to a cell in the domain, + // which has only a node or an edge on the border (not a shared face) + + for (int idomain = 0; idomain < theElems.size(); idomain++) + { + const TIDSortedElemSet& domain = theElems[idomain]; + itface = faceDomains.begin(); + for (; itface != faceDomains.end(); ++itface) + { + std::map 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); + std::set::iterator itn = oldNodes.begin(); + for (; itn != oldNodes.end(); ++itn) + { + int oldId = *itn; + //MESSAGE(" node " << oldId); + std::set cells; + cells.clear(); + vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId); + for (int i=0; iFindElement(GetMeshDS()->fromVtkToSmds(vtkId)); + if (!domain.count(anElem)) + continue; + int vtkType = grid->GetCellType(vtkId); + int downId = grid->CellIdToDownId(vtkId); + if (downId < 0) + { + MESSAGE("doubleNodesOnGroupBoundaries: internal algorithm problem"); + continue; // not OK at this stage of the algorithm: + //no cells created after BuildDownWardConnectivity + } + DownIdType aCell(downId, vtkType); + if (celldom.count(vtkId)) + continue; + cellDomains[aCell][idomain] = vtkId; + celldom[vtkId] = idomain; + } + } + } + } - // --- for each shared face, get the nodes + // --- explore the shared faces domain by domain, to duplicate the nodes in a coherent way + // for each shared face, get the nodes // for each node, for each domain of the face, create a clone of the node - std::map, DownIdCompare>::iterator itface = faceDomains.begin(); - for( ; itface != faceDomains.end();++itface ) + // --- edges at the intersection of 3 or 4 domains, with the order of domains to build + // junction elements of type prism or hexa. the key is the pair of nodesId (lower first) + // the value is the ordered domain ids. (more than 4 domains not taken into account) + + std::map, std::vector > edgesMultiDomains; // nodes of edge --> ordered domains + std::map > mutipleNodes; // nodes muti domains with domain order + + for (int idomain = 0; idomain < theElems.size(); idomain++) { - DownIdType face = itface->first; - std::map domvol = itface->second; - std::set oldNodes; - oldNodes.clear(); - grid->GetNodeIds(oldNodes, face.cellId, face.cellType); - std::set::iterator itn = oldNodes.begin(); - for (;itn != oldNodes.end(); ++itn) + itface = faceDomains.begin(); + for (; itface != faceDomains.end(); ++itface) { - int oldId = *itn; - if (!nodeDomains.count(oldId)) - nodeDomains[oldId] = emptyMap; // create an empty entry for node - std::map::iterator itdom = domvol.begin(); - for(; itdom != domvol.end(); ++itdom) + std::map 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); + bool isMultipleDetected = false; + std::set::iterator itn = oldNodes.begin(); + for (; itn != oldNodes.end(); ++itn) { - int idom = itdom->first; - if ( nodeDomains[oldId].empty() ) - nodeDomains[oldId][idom] = oldId; // keep the old node in the first domain - else + int oldId = *itn; + //MESSAGE(" 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 + std::map::iterator itdom = domvol.begin(); + for (; itdom != domvol.end(); ++itdom) + { + int idom = itdom->first; + //MESSAGE(" domain " << idom); + if (!nodeDomains[oldId].count(idom)) // --- node to clone + { + if (nodeDomains[oldId].size() >= 2) // a multiple node + { + vector orderedDoms; + //MESSAGE("multiple node " << oldId); + isMultipleDetected =true; + if (mutipleNodes.count(oldId)) + orderedDoms = mutipleNodes[oldId]; + else + { + map::iterator it = nodeDomains[oldId].begin(); + for (; it != nodeDomains[oldId].end(); ++it) + orderedDoms.push_back(it->first); + } + orderedDoms.push_back(idom); // TODO order ==> push_front or back + //stringstream txt; + //for (int i=0; iGetPoint(oldId); + SMDS_MeshNode *newNode = meshDS->AddNode(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; + } + } + } + if (isMultipleDetected) // 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; + unsigned char cellType = itface->first.cellType; + int nbEdges = grid->getDownArray(cellType)->getNumberOfDownCells(downId); + const int *downEdgeIds = grid->getDownArray(cellType)->getDownCells(downId); + const unsigned char* edgeType = grid->getDownArray(cellType)->getDownTypes(downId); + for (int ie =0; ie < nbEdges; ie++) { - double *coords = grid->GetPoint(oldId); - SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]); - int newId = newNode->getVtkId(); - nodeDomains[oldId][idom] = newId; // cloned node for other domains + int nodes[3]; + int nbNodes = grid->getDownArray(edgeType[ie])->getNodes(downEdgeIds[ie], nodes); + if (mutipleNodes.count(nodes[0]) && mutipleNodes.count(nodes[nbNodes-1])) + { + vector vn0 = mutipleNodes[nodes[0]]; + vector vn1 = mutipleNodes[nodes[nbNodes - 1]]; + sort( vn0.begin(), vn0.end() ); + sort( vn1.begin(), vn1.end() ); + if (vn0 == vn1) + { + //MESSAGE(" detect edgesMultiDomains " << nodes[0] << " " << nodes[nbNodes - 1]); + double *coords = grid->GetPoint(nodes[0]); + gp_Pnt p0(coords[0], coords[1], coords[2]); + coords = grid->GetPoint(nodes[nbNodes - 1]); + gp_Pnt p1(coords[0], coords[1], coords[2]); + gp_Pnt gref; + int vtkVolIds[1000]; // an edge can belong to a lot of volumes + map 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++) + { + int idom = vn0[id]; + for (int ivol=0; ivolfromVtkToSmds(vtkVolIds[ivol]); + SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId); + if (theElems[idom].count(elem)) + { + SMDS_VtkVolume* svol = dynamic_cast(elem); + domvol[idom] = svol; + //MESSAGE(" domain " << idom << " volume " << elem->GetID()); + double values[3]; + vtkIdType npts = 0; + vtkIdType* pts = 0; + grid->GetCellPoints(vtkVolIds[ivol], npts, pts); + SMDS_VtkVolume::gravityCenter(grid, pts, npts, values); + if (id ==0) + { + gref.SetXYZ(gp_XYZ(values[0], values[1], values[2])); + angleDom[idom] = 0; + } + else + { + gp_Pnt g(values[0], values[1], values[2]); + angleDom[idom] = OrientedAngle(p0, p1, gref, g); // -pisecond << " angle " << ib->first); + } + for (int ino = 0; ino < nbNodes; ino++) + vnodes.push_back(nodes[ino]); + edgesMultiDomains[vnodes] = vdom; // nodes vector --> ordered domains + } + } } } } @@ -10668,42 +11144,128 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector ((domain1 X domain2) --> newId) + // (domain1 X domain2) = domain1 + MAXINT*domain2 + + std::map > nodeQuadDomains; + std::map mapOfJunctionGroups; + if (createJointElems) { itface = faceDomains.begin(); - for( ; itface != faceDomains.end();++itface ) + for (; itface != faceDomains.end(); ++itface) { DownIdType face = itface->first; std::set oldNodes; std::set::iterator itn; oldNodes.clear(); grid->GetNodeIds(oldNodes, face.cellId, face.cellType); - std::map localClonedNodeIds; - std::map domvol = itface->second; - std::map::iterator itdom = domvol.begin(); + std::map domvol = itface->second; + std::map::iterator itdom = domvol.begin(); int dom1 = itdom->first; int vtkVolId = itdom->second; itdom++; int dom2 = itdom->first; + SMDS_MeshVolume *vol = grid->extrudeVolumeFromFace(vtkVolId, dom1, dom2, oldNodes, nodeDomains, + nodeQuadDomains); + stringstream grpname; + grpname << "j_"; + if (dom1 < dom2) + grpname << dom1 << "_" << dom2; + else + grpname << dom2 << "_" << dom1; + int idg; + string namegrp = grpname.str(); + if (!mapOfJunctionGroups.count(namegrp)) + mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg); + SMESHDS_Group *sgrp = dynamic_cast(mapOfJunctionGroups[namegrp]->GetGroupDS()); + if (sgrp) + sgrp->Add(vol->GetID()); + } + } - localClonedNodeIds.clear(); - for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn) + // --- create volumes on multiple domain intersection if requested + // iterate on edgesMultiDomains + + if (createJointElems) + { + std::map, std::vector >::iterator ite = edgesMultiDomains.begin(); + for (; ite != edgesMultiDomains.end(); ++ite) + { + vector nodes = ite->first; + vector orderDom = ite->second; + vector orderedNodes; + if (nodes.size() == 2) { - int oldId = *itn; - int refid = oldId; - if (nodeDomains[oldId].count(dom1)) - refid = nodeDomains[oldId][dom1]; - else - MESSAGE("--- problem domain node " << dom1 << " " << oldId); - int newid = oldId; - if (nodeDomains[oldId].count(dom2)) - newid = nodeDomains[oldId][dom2]; - else - MESSAGE("--- problem domain node " << dom2 << " " << oldId); - localClonedNodeIds[oldId] = newid; + //MESSAGE(" use edgesMultiDomains " << nodes[0] << " " << nodes[1]); + for (int ino=0; ino < nodes.size(); ino++) + if (orderDom.size() == 3) + for (int idom = 0; idom =0; idom--) + orderedNodes.push_back( nodeDomains[nodes[ino]][orderDom[idom]] ); + SMDS_MeshVolume* vol = this->GetMeshDS()->AddVolumeFromVtkIds(orderedNodes); + + stringstream grpname; + grpname << "mj_"; + grpname << 0 << "_" << 0; + int idg; + string namegrp = grpname.str(); + if (!mapOfJunctionGroups.count(namegrp)) + mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg); + SMESHDS_Group *sgrp = dynamic_cast(mapOfJunctionGroups[namegrp]->GetGroupDS()); + if (sgrp) + sgrp->Add(vol->GetID()); + } + else + { + //MESSAGE("Quadratic multiple joints not implemented"); + // TODO quadratic nodes + } + } + } + + // --- list the explicit faces and edges of the mesh that need to be modified, + // i.e. faces and edges built with one or more duplicated nodes. + // associate these faces or edges to their corresponding domain. + // only the first domain found is kept when a face or edge is shared + + std::map, DownIdCompare> faceOrEdgeDom; // cellToModify --> (id domain --> id cell) + std::map feDom; // vtk id of cell to modify --> id domain + faceOrEdgeDom.clear(); + feDom.clear(); + + for (int idomain = 0; idomain < theElems.size(); idomain++) + { + std::map >::const_iterator itnod = nodeDomains.begin(); + for (; itnod != nodeDomains.end(); ++itnod) + { + int oldId = itnod->first; + //MESSAGE(" node " << oldId); + vtkCellLinks::Link l = grid->GetCellLinks()->GetLink(oldId); + for (int i = 0; i < l.ncells; i++) + { + int vtkId = l.cells[i]; + int vtkType = grid->GetCellType(vtkId); + int downId = grid->CellIdToDownId(vtkId); + if (downId < 0) + continue; // new cells: not to be modified + DownIdType aCell(downId, vtkType); + int volParents[1000]; + int nbvol = grid->GetParentVolumes(volParents, vtkId); + for (int j = 0; j < nbvol; j++) + if (celldom.count(volParents[j]) && (celldom[volParents[j]] == idomain)) + if (!feDom.count(vtkId)) + { + feDom[vtkId] = idomain; + faceOrEdgeDom[aCell] = emptyMap; + faceOrEdgeDom[aCell][idomain] = vtkId; // affect face or edge to the first domain only + //MESSAGE("affect cell " << this->GetMeshDS()->fromVtkToSmds(vtkId) << " domain " << idomain + // << " type " << vtkType << " downId " << downId); + } } - meshDS->extrudeVolumeFromFace(vtkVolId, localClonedNodeIds); } } @@ -10711,40 +11273,212 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector, DownIdCompare>* maps[3] = {&faceDomains, &cellDomains, &faceOrEdgeDom}; + for (int m=0; m<3; m++) { - DownIdType face = itface->first; - std::set oldNodes; - std::set::iterator itn; - oldNodes.clear(); - grid->GetNodeIds(oldNodes, face.cellId, face.cellType); - std::map localClonedNodeIds; - - std::map domvol = itface->second; - std::map::iterator itdom = domvol.begin(); - for(; itdom != domvol.end(); ++itdom) + std::map, DownIdCompare>* amap = maps[m]; + itface = (*amap).begin(); + for (; itface != (*amap).end(); ++itface) { - int idom = itdom->first; - int vtkVolId = itdom->second; - localClonedNodeIds.clear(); - for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn) + DownIdType face = itface->first; + std::set oldNodes; + std::set::iterator itn; + oldNodes.clear(); + grid->GetNodeIds(oldNodes, face.cellId, face.cellType); + //MESSAGE("examine cell, downId " << face.cellId << " type " << int(face.cellType)); + std::map localClonedNodeIds; + + std::map domvol = itface->second; + std::map::iterator itdom = domvol.begin(); + for (; itdom != domvol.end(); ++itdom) { - int oldId = *itn; - if (nodeDomains[oldId].count(idom)) - localClonedNodeIds[oldId] = nodeDomains[oldId][idom]; + int idom = itdom->first; + int vtkVolId = itdom->second; + //MESSAGE("modify nodes of cell " << this->GetMeshDS()->fromVtkToSmds(vtkVolId) << " domain " << idom); + localClonedNodeIds.clear(); + for (itn = oldNodes.begin(); itn != oldNodes.end(); ++itn) + { + int oldId = *itn; + if (nodeDomains[oldId].count(idom)) + { + localClonedNodeIds[oldId] = nodeDomains[oldId][idom]; + //MESSAGE(" node " << oldId << " --> " << localClonedNodeIds[oldId]); + } + } + meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds); } - meshDS->ModifyCellNodes(vtkVolId, localClonedNodeIds); } } + + meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory grid->BuildLinks(); - // TODO replace also old nodes by new nodes in faces and edges CHRONOSTOP(50); counters::stats(); return true; } +/*! + * \brief Double nodes on some external faces and create flat elements. + * Flat elements are mainly used by some types of mechanic calculations. + * + * Each group of the list must be constituted of faces. + * Triangles are transformed in prisms, and quadrangles in hexahedrons. + * @param theElems - list of groups of faces, where a group of faces is a set of + * SMDS_MeshElements sorted by Id. + * @return TRUE if operation has been completed successfully, FALSE otherwise + */ +bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vector& theElems) +{ + MESSAGE("-------------------------------------------------"); + MESSAGE("SMESH_MeshEditor::CreateFlatElementsOnFacesGroups"); + MESSAGE("-------------------------------------------------"); + + SMESHDS_Mesh *meshDS = this->myMesh->GetMeshDS(); + + // --- For each group of faces + // duplicate the nodes, create a flat element based on the face + // replace the nodes of the faces by their clones + + std::map clonedNodes; + std::map intermediateNodes; + clonedNodes.clear(); + intermediateNodes.clear(); + std::map mapOfJunctionGroups; + mapOfJunctionGroups.clear(); + + for (int idom = 0; idom < theElems.size(); idom++) + { + const TIDSortedElemSet& domain = theElems[idom]; + TIDSortedElemSet::const_iterator elemItr = domain.begin(); + for (; elemItr != domain.end(); ++elemItr) + { + SMDS_MeshElement* anElem = (SMDS_MeshElement*) *elemItr; + SMDS_MeshFace* aFace = dynamic_cast (anElem); + if (!aFace) + continue; + // MESSAGE("aFace=" << aFace->GetID()); + bool isQuad = aFace->IsQuadratic(); + vector ln0, ln1, ln2, ln3, ln4; + + // --- clone the nodes, create intermediate nodes for non medium nodes of a quad face + + SMDS_ElemIteratorPtr nodeIt = aFace->nodesIterator(); + while (nodeIt->more()) + { + const SMDS_MeshNode* node = static_cast (nodeIt->next()); + bool isMedium = isQuad && (aFace->IsMediumNode(node)); + if (isMedium) + ln2.push_back(node); + else + ln0.push_back(node); + + const SMDS_MeshNode* clone = 0; + if (!clonedNodes.count(node)) + { + clone = meshDS->AddNode(node->X(), node->Y(), node->Z()); + clonedNodes[node] = clone; + } + else + clone = clonedNodes[node]; + + if (isMedium) + ln3.push_back(clone); + else + ln1.push_back(clone); + + const SMDS_MeshNode* inter = 0; + if (isQuad && (!isMedium)) + { + if (!intermediateNodes.count(node)) + { + inter = meshDS->AddNode(node->X(), node->Y(), node->Z()); + intermediateNodes[node] = inter; + } + else + inter = intermediateNodes[node]; + ln4.push_back(inter); + } + } + + // --- extrude the face + + vector ln; + SMDS_MeshVolume* vol = 0; + vtkIdType aType = aFace->GetVtkType(); + switch (aType) + { + case VTK_TRIANGLE: + vol = meshDS->AddVolume(ln0[2], ln0[1], ln0[0], ln1[2], ln1[1], ln1[0]); + // MESSAGE("vol prism " << vol->GetID()); + ln.push_back(ln1[0]); + ln.push_back(ln1[1]); + ln.push_back(ln1[2]); + break; + case VTK_QUAD: + vol = meshDS->AddVolume(ln0[3], ln0[2], ln0[1], ln0[0], ln1[3], ln1[2], ln1[1], ln1[0]); + // MESSAGE("vol hexa " << vol->GetID()); + ln.push_back(ln1[0]); + ln.push_back(ln1[1]); + ln.push_back(ln1[2]); + ln.push_back(ln1[3]); + break; + case VTK_QUADRATIC_TRIANGLE: + vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln0[0], ln0[1], ln0[2], ln3[0], ln3[1], ln3[2], + ln2[0], ln2[1], ln2[2], ln4[0], ln4[1], ln4[2]); + // MESSAGE("vol quad prism " << vol->GetID()); + ln.push_back(ln1[0]); + ln.push_back(ln1[1]); + ln.push_back(ln1[2]); + ln.push_back(ln3[0]); + ln.push_back(ln3[1]); + ln.push_back(ln3[2]); + break; + case VTK_QUADRATIC_QUAD: +// vol = meshDS->AddVolume(ln0[0], ln0[1], ln0[2], ln0[3], ln1[0], ln1[1], ln1[2], ln1[3], +// ln2[0], ln2[1], ln2[2], ln2[3], ln3[0], ln3[1], ln3[2], ln3[3], +// ln4[0], ln4[1], ln4[2], ln4[3]); + vol = meshDS->AddVolume(ln1[0], ln1[1], ln1[2], ln1[3], ln0[0], ln0[1], ln0[2], ln0[3], + ln3[0], ln3[1], ln3[2], ln3[3], ln2[0], ln2[1], ln2[2], ln2[3], + ln4[0], ln4[1], ln4[2], ln4[3]); + // MESSAGE("vol quad hexa " << vol->GetID()); + ln.push_back(ln1[0]); + ln.push_back(ln1[1]); + ln.push_back(ln1[2]); + ln.push_back(ln1[3]); + ln.push_back(ln3[0]); + ln.push_back(ln3[1]); + ln.push_back(ln3[2]); + ln.push_back(ln3[3]); + break; + case VTK_POLYGON: + break; + default: + break; + } + + if (vol) + { + stringstream grpname; + grpname << "jf_"; + grpname << idom; + int idg; + string namegrp = grpname.str(); + if (!mapOfJunctionGroups.count(namegrp)) + mapOfJunctionGroups[namegrp] = this->myMesh->AddGroup(SMDSAbs_Volume, namegrp.c_str(), idg); + SMESHDS_Group *sgrp = dynamic_cast(mapOfJunctionGroups[namegrp]->GetGroupDS()); + if (sgrp) + sgrp->Add(vol->GetID()); + } + + // --- modify the face + + aFace->ChangeNodes(&ln[0], ln.size()); + } + } + return true; +} + //================================================================================ /*! * \brief Generates skin mesh (containing 2D cells) from 3D mesh @@ -10815,15 +11549,23 @@ namespace * \param toCopyElements - if true, the checked elements will be copied into the targetMesh * \param toCopyExistingBondary - if true, not only new but also pre-existing * boundary elements will be copied into the targetMesh + * \param toAddExistingBondary - if true, not only new but also pre-existing + * boundary elements will be added into the new group + * \param aroundElements - if true, elements will be created on boundary of given + * elements else, on boundary of the whole mesh. This + * option works for 2D elements only. + * \return nb of added boundary elements */ //================================================================================ -void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, - Bnd_Dimension dimension, - SMESH_Group* group/*=0*/, - SMESH_Mesh* targetMesh/*=0*/, - bool toCopyElements/*=false*/, - bool toCopyExistingBondary/*=false*/) +int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, + Bnd_Dimension dimension, + SMESH_Group* group/*=0*/, + SMESH_Mesh* targetMesh/*=0*/, + bool toCopyElements/*=false*/, + bool toCopyExistingBondary/*=false*/, + bool toAddExistingBondary/*= false*/, + bool aroundElements/*= false*/) { SMDSAbs_ElementType missType = (dimension == BND_2DFROM3D) ? SMDSAbs_Face : SMDSAbs_Edge; SMDSAbs_ElementType elemType = (dimension == BND_1DFROM2D) ? SMDSAbs_Face : SMDSAbs_Volume; @@ -10831,14 +11573,24 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, if ( !elements.empty() && (*elements.begin())->GetType() != elemType ) throw SALOME_Exception(LOCALIZED("wrong element type")); + if ( aroundElements && elemType == SMDSAbs_Volume ) + throw SALOME_Exception(LOCALIZED("wrong element type for aroundElements==true")); + if ( !targetMesh ) toCopyElements = toCopyExistingBondary = false; SMESH_MeshEditor tgtEditor( targetMesh ? targetMesh : myMesh ); SMESHDS_Mesh* aMesh = GetMeshDS(), *tgtMeshDS = tgtEditor.GetMeshDS(); + int nbAddedBnd = 0; + + // editor adding present bnd elements and optionally holding elements to add to the group + SMESH_MeshEditor* presentEditor; + SMESH_MeshEditor tgtEditor2( tgtEditor.GetMesh() ); + presentEditor = toAddExistingBondary ? &tgtEditor : &tgtEditor2; SMDS_VolumeTool vTool; - TIDSortedElemSet emptySet, avoidSet; + TIDSortedElemSet avoidSet; + const TIDSortedElemSet emptySet, *elemSet = aroundElements ? &elements : &emptySet; int inode; typedef vector TConnectivity; @@ -10854,7 +11606,9 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, const SMDS_MeshElement* elem = eIt->next(); 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; @@ -10906,7 +11660,7 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, { nodes[0] = elem->GetNode(i); nodes[1] = elem->GetNode((i+1)%nbNodes); - if ( FindFaceInSet( nodes[0], nodes[1], emptySet, avoidSet)) + if ( FindFaceInSet( nodes[0], nodes[1], *elemSet, avoidSet)) continue; // not free link //if ( iQuad ) @@ -10919,7 +11673,9 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, } } + // --------------------------------- // 2. Add missing boundary elements + // --------------------------------- if ( targetMesh != myMesh ) // instead of making a map of nodes in this mesh and targetMesh, // we create nodes with same IDs. We can renumber them later, if needed @@ -10929,16 +11685,28 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, TConnectivity nodes( srcNodes.size() ); for ( inode = 0; inode < nodes.size(); ++inode ) nodes[inode] = getNodeWithSameID( tgtMeshDS, srcNodes[inode] ); + if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes, + missType, + /*noMedium=*/true)) + continue; tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4); + ++nbAddedBnd; } else for ( int i = 0; i < missingBndElems.size(); ++i ) { - TConnectivity& nodes = missingBndElems[i]; + TConnectivity& nodes = missingBndElems[i]; + if ( aroundElements && tgtEditor.GetMeshDS()->FindElement( nodes, + missType, + /*noMedium=*/true)) + continue; tgtEditor.AddElement(nodes, missType, elem->IsPoly() && nodes.size()/(iQuad+1)>4); + ++nbAddedBnd; } + // ---------------------------------- // 3. Copy present boundary elements + // ---------------------------------- if ( toCopyExistingBondary ) for ( int i = 0 ; i < presentBndElems.size(); ++i ) { @@ -10946,13 +11714,19 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, TConnectivity nodes( e->NbNodes() ); for ( inode = 0; inode < nodes.size(); ++inode ) nodes[inode] = getNodeWithSameID( tgtMeshDS, e->GetNode(inode) ); - tgtEditor.AddElement(nodes, missType, e->IsPoly()); - // leave only missing elements in tgtEditor.myLastCreatedElems - tgtEditor.myLastCreatedElems.Remove( tgtEditor.myLastCreatedElems.Size() ); + presentEditor->AddElement(nodes, missType, e->IsPoly()); + } + else // store present elements to add them to a group + for ( int i = 0 ; i < presentBndElems.size(); ++i ) + { + presentEditor->myLastCreatedElems.Append(presentBndElems[i]); } + } // loop on given elements - // 4. Fill group with missing boundary elements + // --------------------------------------------- + // 4. Fill group with boundary elements + // --------------------------------------------- if ( group ) { if ( SMESHDS_Group* g = dynamic_cast( group->GetGroupDS() )) @@ -10960,9 +11734,12 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, g->SMDSGroup().Add( tgtEditor.myLastCreatedElems( i+1 )); } tgtEditor.myLastCreatedElems.Clear(); + tgtEditor2.myLastCreatedElems.Clear(); + // ----------------------- // 5. Copy given elements - if ( toCopyElements ) + // ----------------------- + if ( toCopyElements && targetMesh != myMesh ) { if (elements.empty()) eIt = aMesh->elementsIterator(elemType); @@ -10979,5 +11756,5 @@ void SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, tgtEditor.myLastCreatedElems.Clear(); } } - return; + return nbAddedBnd; }