X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MeshEditor.cxx;h=7083261a9226c5625629d12ffd14c0b1a06e74b5;hp=49d3e14c7ed6f452ce8e775e72c38f5c489806fb;hb=2e9f6a1d3399b1ea9b366f969e81c725a5a5a628;hpb=094287b4dfbd30d6b17cb9d6cb892e62388c2170 diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index 49d3e14c7..7083261a9 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -6,7 +6,7 @@ // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // License as published by the Free Software Foundation; either -// version 2.1 of the License. +// version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of @@ -50,6 +50,7 @@ #include #include "utilities.h" +#include "chrono.hxx" #include #include @@ -518,14 +519,12 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem) } 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; + SMESHDS_SubMeshIteratorPtr smIt = GetMeshDS()->SubMeshes(); + while ( const SMESHDS_SubMesh* sm = smIt->next() ) + if ( sm->Contains( theElem )) + return sm->GetID(); } - //MESSAGE ("::FindShape() - SHAPE NOT FOUND") return 0; } @@ -658,7 +657,7 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1, // put nodes in array and find out indices of the same ones const SMDS_MeshNode* aNodes [6]; - int sameInd [] = { 0, 0, 0, 0, 0, 0 }; + int sameInd [] = { -1, -1, -1, -1, -1, -1 }; int i = 0; SMDS_ElemIteratorPtr it = theTria1->nodesIterator(); while ( it->more() ) { @@ -684,15 +683,15 @@ bool SMESH_MeshEditor::InverseDiag (const SMDS_MeshElement * theTria1, } // find indices of 1,2 and of A,B in theTria1 - int iA = 0, iB = 0, i1 = 0, i2 = 0; + int iA = -1, iB = 0, i1 = 0, i2 = 0; for ( i = 0; i < 6; i++ ) { - if ( sameInd [ i ] == 0 ) { + if ( sameInd [ i ] == -1 ) { if ( i < 3 ) i1 = i; else i2 = i; } else if (i < 3) { - if ( iA ) iB = i; - else iA = i; + if ( iA >= 0) iB = i; + else iA = i; } } // nodes 1 and 2 should not be the same @@ -1174,7 +1173,7 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, avoidSet.clear(); avoidSet.insert(theFace); - NLink link( theFace->GetNode( 0 ), 0 ); + NLink link( theFace->GetNode( 0 ), (SMDS_MeshNode *) 0 ); const int nbNodes = theFace->NbCornerNodes(); for ( int i = 0; i < nbNodes; ++i ) // loop on links of theFace @@ -1245,6 +1244,92 @@ int SMESH_MeshEditor::Reorient2D (TIDSortedElemSet & theFaces, return nbReori; } +//================================================================================ +/*! + * \brief Reorient faces basing on orientation of adjacent volumes. + * \param theFaces - faces to reorient. If empty, all mesh faces are treated. + * \param theVolumes - reference volumes. + * \param theOutsideNormal - to orient faces to have their normal + * pointing either \a outside or \a inside the adjacent volumes. + * \return number of reoriented faces. + */ +//================================================================================ + +int SMESH_MeshEditor::Reorient2DBy3D (TIDSortedElemSet & theFaces, + TIDSortedElemSet & theVolumes, + const bool theOutsideNormal) +{ + int nbReori = 0; + + SMDS_ElemIteratorPtr faceIt; + if ( theFaces.empty() ) + faceIt = GetMeshDS()->elementsIterator( SMDSAbs_Face ); + else + faceIt = elemSetIterator( theFaces ); + + vector< const SMDS_MeshNode* > faceNodes; + TIDSortedElemSet checkedVolumes; + set< const SMDS_MeshNode* > faceNodesSet; + SMDS_VolumeTool volumeTool; + + while ( faceIt->more() ) // loop on given faces + { + const SMDS_MeshElement* face = faceIt->next(); + if ( face->GetType() != SMDSAbs_Face ) + continue; + + const int nbCornersNodes = face->NbCornerNodes(); + faceNodes.assign( face->begin_nodes(), face->end_nodes() ); + + checkedVolumes.clear(); + SMDS_ElemIteratorPtr vIt = faceNodes[ 0 ]->GetInverseElementIterator( SMDSAbs_Volume ); + while ( vIt->more() ) + { + const SMDS_MeshElement* volume = vIt->next(); + + if ( !checkedVolumes.insert( volume ).second ) + continue; + if ( !theVolumes.empty() && !theVolumes.count( volume )) + continue; + + // is volume adjacent? + bool allNodesCommon = true; + for ( int iN = 1; iN < nbCornersNodes && allNodesCommon; ++iN ) + allNodesCommon = ( volume->GetNodeIndex( faceNodes[ iN ]) > -1 ); + if ( !allNodesCommon ) + continue; + + // get nodes of a corresponding volume facet + faceNodesSet.clear(); + faceNodesSet.insert( faceNodes.begin(), faceNodes.end() ); + volumeTool.Set( volume ); + int facetID = volumeTool.GetFaceIndex( faceNodesSet ); + if ( facetID < 0 ) continue; + volumeTool.SetExternalNormal(); + const SMDS_MeshNode** facetNodes = volumeTool.GetFaceNodes( facetID ); + + // compare order of faceNodes and facetNodes + const int iQ = 1 + ( nbCornersNodes < faceNodes.size() ); + int iNN[2]; + for ( int i = 0; i < 2; ++i ) + { + const SMDS_MeshNode* n = facetNodes[ i*iQ ]; + for ( int iN = 0; iN < nbCornersNodes; ++iN ) + if ( faceNodes[ iN ] == n ) + { + iNN[ i ] = iN; + break; + } + } + bool isOutside = Abs( iNN[0]-iNN[1] ) == 1 ? iNN[0] < iNN[1] : iNN[0] > iNN[1]; + if ( isOutside != theOutsideNormal ) + nbReori += Reorient( face ); + } + } // loop on given faces + + return nbReori; +} + //======================================================================= //function : getBadRate //purpose : @@ -1603,44 +1688,110 @@ namespace const int* thePentaTo3[6] = { thePentaTo3_1, thePentaTo3_2, thePentaTo3_3, thePentaTo3_4, thePentaTo3_5, thePentaTo3_6 }; + // Methods of splitting hexahedron into prisms + + const int theHexTo4Prisms_BT[6*4+1] = // bottom-top + { + 0, 1, 8, 4, 5, 9, 1, 2, 8, 5, 6, 9, 2, 3, 8, 6, 7, 9, 3, 0, 8, 7, 4, 9, -1 + }; + const int theHexTo4Prisms_LR[6*4+1] = // left-right + { + 1, 0, 8, 2, 3, 9, 0, 4, 8, 3, 7, 9, 4, 5, 8, 7, 6, 9, 5, 1, 8, 6, 2, 9, -1 + }; + const int theHexTo4Prisms_FB[6*4+1] = // front-back + { + 0, 3, 9, 1, 2, 8, 3, 7, 9, 2, 6, 8, 7, 4, 9, 6, 5, 8, 4, 0, 9, 5, 1, 8, -1 + }; + + const int theHexTo2Prisms_BT_1[6*2+1] = + { + 0, 1, 3, 4, 5, 7, 1, 2, 3, 5, 6, 7, -1 + }; + const int theHexTo2Prisms_BT_2[6*2+1] = + { + 0, 1, 2, 4, 5, 6, 0, 2, 3, 4, 6, 7, -1 + }; + const int* theHexTo2Prisms_BT[2] = { theHexTo2Prisms_BT_1, theHexTo2Prisms_BT_2 }; + + const int theHexTo2Prisms_LR_1[6*2+1] = + { + 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1 + }; + const int theHexTo2Prisms_LR_2[6*2+1] = + { + 1, 0, 4, 2, 3, 7, 1, 4, 5, 2, 7, 6, -1 + }; + const int* theHexTo2Prisms_LR[2] = { theHexTo2Prisms_LR_1, theHexTo2Prisms_LR_2 }; + + const int theHexTo2Prisms_FB_1[6*2+1] = + { + 0, 3, 4, 1, 2, 5, 3, 7, 4, 2, 6, 5, -1 + }; + const int theHexTo2Prisms_FB_2[6*2+1] = + { + 0, 3, 7, 1, 2, 7, 0, 7, 4, 1, 6, 5, -1 + }; + const int* theHexTo2Prisms_FB[2] = { theHexTo2Prisms_FB_1, theHexTo2Prisms_FB_2 }; + + struct TTriangleFacet //!< stores indices of three nodes of tetra facet { int _n1, _n2, _n3; TTriangleFacet(int n1, int n2, int n3): _n1(n1), _n2(n2), _n3(n3) {} bool contains(int n) const { return ( n == _n1 || n == _n2 || n == _n3 ); } - bool hasAdjacentTetra( const SMDS_MeshElement* elem ) const; + bool hasAdjacentVol( const SMDS_MeshElement* elem, + const SMDSAbs_GeometryType geom = SMDSGeom_TETRA) const; }; struct TSplitMethod { - int _nbTetra; + int _nbSplits; + int _nbCorners; const int* _connectivity; //!< foursomes of tetra connectivy finished by -1 bool _baryNode; //!< additional node is to be created at cell barycenter bool _ownConn; //!< to delete _connectivity in destructor map _faceBaryNode; //!< map face index to node at BC of face TSplitMethod( int nbTet=0, const int* conn=0, bool addNode=false) - : _nbTetra(nbTet), _connectivity(conn), _baryNode(addNode), _ownConn(false) {} + : _nbSplits(nbTet), _nbCorners(4), _connectivity(conn), _baryNode(addNode), _ownConn(false) {} ~TSplitMethod() { if ( _ownConn ) delete [] _connectivity; _connectivity = 0; } bool hasFacet( const TTriangleFacet& facet ) const { - const int* tetConn = _connectivity; - for ( ; tetConn[0] >= 0; tetConn += 4 ) - if (( facet.contains( tetConn[0] ) + - facet.contains( tetConn[1] ) + - facet.contains( tetConn[2] ) + - facet.contains( tetConn[3] )) == 3 ) - return true; + if ( _nbCorners == 4 ) + { + const int* tetConn = _connectivity; + for ( ; tetConn[0] >= 0; tetConn += 4 ) + if (( facet.contains( tetConn[0] ) + + facet.contains( tetConn[1] ) + + facet.contains( tetConn[2] ) + + facet.contains( tetConn[3] )) == 3 ) + return true; + } + else // prism, _nbCorners == 6 + { + const int* prismConn = _connectivity; + for ( ; prismConn[0] >= 0; prismConn += 6 ) + { + if (( facet.contains( prismConn[0] ) && + facet.contains( prismConn[1] ) && + facet.contains( prismConn[2] )) + || + ( facet.contains( prismConn[3] ) && + facet.contains( prismConn[4] ) && + facet.contains( prismConn[5] ))) + return true; + } + } return false; } }; //======================================================================= /*! - * \brief return TSplitMethod for the given element + * \brief return TSplitMethod for the given element to split into tetrahedra */ //======================================================================= - TSplitMethod getSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags) + TSplitMethod getTetraSplitMethod( SMDS_VolumeTool& vol, const int theMethodFlags) { const int iQ = vol.Element()->IsQuadratic() ? 2 : 1; @@ -1665,8 +1816,8 @@ namespace { TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] ); TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] ); - if ( t012.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t012 ); - else if ( t123.hasAdjacentTetra( vol.Element() )) triaSplits.push_back( t123 ); + if ( t012.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t012 ); + else if ( t123.hasAdjacentVol( vol.Element() )) triaSplits.push_back( t123 ); } else { @@ -1679,7 +1830,7 @@ namespace TTriangleFacet t023( nInd[ iQ * ( iCom )], nInd[ iQ * ( (iCom+2)%nbNodes )], nInd[ iQ * ( (iCom+3)%nbNodes )]); - if ( t012.hasAdjacentTetra( vol.Element() ) && t023.hasAdjacentTetra( vol.Element() )) + if ( t012.hasAdjacentVol( vol.Element() ) && t023.hasAdjacentVol( vol.Element() )) { triaSplits.push_back( t012 ); triaSplits.push_back( t023 ); @@ -1719,12 +1870,12 @@ namespace default: nbVariants = 0; } - for ( int variant = 0; variant < nbVariants && method._nbTetra == 0; ++variant ) + for ( int variant = 0; variant < nbVariants && method._nbSplits == 0; ++variant ) { // check method compliancy with adjacent tetras, // all found splits must be among facets of tetras described by this method method = TSplitMethod( nbTet, connVariants[variant] ); - if ( hasAdjacentSplits && method._nbTetra > 0 ) + if ( hasAdjacentSplits && method._nbSplits > 0 ) { bool facetCreated = true; for ( int iF = 0; facetCreated && iF < triaSplitsByFace.size(); ++iF ) @@ -1738,7 +1889,7 @@ namespace } } } - if ( method._nbTetra < 1 ) + if ( method._nbSplits < 1 ) { // No standard method is applicable, use a generic solution: // each facet of a volume is split into triangles and @@ -1832,7 +1983,7 @@ namespace connectivity[ connSize++ ] = baryCenInd; } } - method._nbTetra += nbTet; + method._nbSplits += nbTet; } // loop on volume faces @@ -1842,13 +1993,132 @@ namespace return method; } + //======================================================================= + /*! + * \brief return TSplitMethod to split haxhedron into prisms + */ + //======================================================================= + + TSplitMethod getPrismSplitMethod( SMDS_VolumeTool& vol, + const int methodFlags, + const int facetToSplit) + { + // order of facets in HEX according to SMDS_VolumeTool::Hexa_F : + // B, T, L, B, R, F + const int iF = ( facetToSplit < 2 ) ? 0 : 1 + ( facetToSplit-2 ) % 2; // [0,1,2] + + if ( methodFlags == SMESH_MeshEditor::HEXA_TO_4_PRISMS ) + { + static TSplitMethod to4methods[4]; // order BT, LR, FB + if ( to4methods[iF]._nbSplits == 0 ) + { + switch ( iF ) { + case 0: + to4methods[iF]._connectivity = theHexTo4Prisms_BT; + to4methods[iF]._faceBaryNode[ 0 ] = 0; + to4methods[iF]._faceBaryNode[ 1 ] = 0; + break; + case 1: + to4methods[iF]._connectivity = theHexTo4Prisms_LR; + to4methods[iF]._faceBaryNode[ 2 ] = 0; + to4methods[iF]._faceBaryNode[ 4 ] = 0; + break; + case 2: + to4methods[iF]._connectivity = theHexTo4Prisms_FB; + to4methods[iF]._faceBaryNode[ 3 ] = 0; + to4methods[iF]._faceBaryNode[ 5 ] = 0; + break; + default: return to4methods[3]; + } + to4methods[iF]._nbSplits = 4; + to4methods[iF]._nbCorners = 6; + } + return to4methods[iF]; + } + // else if ( methodFlags == HEXA_TO_2_PRISMS ) + + TSplitMethod method; + + const int iQ = vol.Element()->IsQuadratic() ? 2 : 1; + + const int nbVariants = 2, nbSplits = 2; + const int** connVariants = 0; + switch ( iF ) { + case 0: connVariants = theHexTo2Prisms_BT; break; + case 1: connVariants = theHexTo2Prisms_LR; break; + case 2: connVariants = theHexTo2Prisms_FB; break; + default: return method; + } + + // look for prisms adjacent via facetToSplit and an opposite one + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit; + int nbNodes = vol.NbFaceNodes( iFacet ) / iQ; + if ( nbNodes != 4 ) return method; + + const int* nInd = vol.GetFaceNodesIndices( iFacet ); + TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] ); + TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] ); + TTriangleFacet* t; + if ( t012.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA )) + t = &t012; + else if ( t123.hasAdjacentVol( vol.Element(), SMDSGeom_PENTA )) + t = &t123; + else + continue; + + // there are adjacent prism + for ( int variant = 0; variant < nbVariants; ++variant ) + { + // check method compliancy with adjacent prisms, + // the found prism facets must be among facets of prisms described by current method + method._nbSplits = nbSplits; + method._nbCorners = 6; + method._connectivity = connVariants[ variant ]; + if ( method.hasFacet( *t )) + return method; + } + } + + // No adjacent prisms. Select a variant with a best aspect ratio. + + double badness[2] = { 0, 0 }; + static SMESH::Controls::NumericalFunctorPtr aspectRatio( new SMESH::Controls::AspectRatio); + const SMDS_MeshNode** nodes = vol.GetNodes(); + for ( int variant = 0; variant < nbVariants; ++variant ) + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + int iFacet = is2nd ? vol.GetOppFaceIndexOfHex( facetToSplit ) : facetToSplit; + const int* nInd = vol.GetFaceNodesIndices( iFacet ); + + method._connectivity = connVariants[ variant ]; + TTriangleFacet t012( nInd[0*iQ], nInd[1*iQ], nInd[2*iQ] ); + TTriangleFacet t123( nInd[1*iQ], nInd[2*iQ], nInd[3*iQ] ); + TTriangleFacet* t = ( method.hasFacet( t012 )) ? & t012 : & t123; + + SMDS_FaceOfNodes tria ( nodes[ t->_n1 ], + nodes[ t->_n2 ], + nodes[ t->_n3 ] ); + badness[ variant ] += getBadRate( &tria, aspectRatio ); + } + const int iBetter = ( badness[1] < badness[0] && badness[0]-badness[1] > 0.1 * badness[0] ); + + method._nbSplits = nbSplits; + method._nbCorners = 6; + method._connectivity = connVariants[ iBetter ]; + + return method; + } + //================================================================================ /*! * \brief Check if there is a tetraherdon adjacent to the given element via this facet */ //================================================================================ - bool TTriangleFacet::hasAdjacentTetra( const SMDS_MeshElement* elem ) const + bool TTriangleFacet::hasAdjacentVol( const SMDS_MeshElement* elem, + const SMDSAbs_GeometryType geom ) const { // find the tetrahedron including the three nodes of facet const SMDS_MeshNode* n1 = elem->GetNode(_n1); @@ -1858,16 +2128,16 @@ namespace while ( volIt1->more() ) { const SMDS_MeshElement* v = volIt1->next(); - SMDSAbs_EntityType type = v->GetEntityType(); - if ( type != SMDSEntity_Tetra && type != SMDSEntity_Quad_Tetra ) + if ( v->GetGeomType() != geom ) continue; - if ( type == SMDSEntity_Quad_Tetra && v->GetNodeIndex( n1 ) > 3 ) + const int lastCornerInd = v->NbCornerNodes() - 1; + if ( v->IsQuadratic() && v->GetNodeIndex( n1 ) > lastCornerInd ) continue; // medium node not allowed const int ind2 = v->GetNodeIndex( n2 ); - if ( ind2 < 0 || 3 < ind2 ) + if ( ind2 < 0 || lastCornerInd < ind2 ) continue; const int ind3 = v->GetNodeIndex( n3 ); - if ( ind3 < 0 || 3 < ind3 ) + if ( ind3 < 0 || lastCornerInd < ind3 ) continue; return true; } @@ -1900,19 +2170,23 @@ namespace } // namespace //======================================================================= -//function : SplitVolumesIntoTetra -//purpose : Split volume elements into tetrahedra. +//function : SplitVolumes +//purpose : Split volume elements into tetrahedra or prisms. +// If facet ID < 0, element is split into tetrahedra, +// else a hexahedron is split into prisms so that the given facet is +// split into triangles //======================================================================= -void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, - const int theMethodFlags) +void SMESH_MeshEditor::SplitVolumes (const TFacetOfElem & theElems, + const int theMethodFlags) { // std-like iterator on coordinates of nodes of mesh element typedef SMDS_StdIterator< SMESH_TNodeXYZ, SMDS_ElemIteratorPtr > NXyzIterator; NXyzIterator xyzEnd; SMDS_VolumeTool volTool; - SMESH_MesherHelper helper( *GetMesh()); + SMESH_MesherHelper helper( *GetMesh()), fHelper(*GetMesh()); + fHelper.ToFixNodeParameters( true ); SMESHDS_SubMesh* subMesh = 0;//GetMeshDS()->MeshElements(1); SMESHDS_SubMesh* fSubMesh = 0;//subMesh; @@ -1923,29 +2197,33 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, map< TVolumeFaceKey, const SMDS_MeshNode* > volFace2BaryNode; double bc[3]; - TIDSortedElemSet::const_iterator elem = theElems.begin(); - for ( ; elem != theElems.end(); ++elem ) + TFacetOfElem::const_iterator elem2facet = theElems.begin(); + for ( ; elem2facet != theElems.end(); ++elem2facet ) { - if ( (*elem)->GetType() != SMDSAbs_Volume ) + const SMDS_MeshElement* elem = elem2facet->first; + const int facetToSplit = elem2facet->second; + if ( elem->GetType() != SMDSAbs_Volume ) continue; - SMDSAbs_EntityType geomType = (*elem)->GetEntityType(); + const SMDSAbs_EntityType geomType = elem->GetEntityType(); if ( geomType == SMDSEntity_Tetra || geomType == SMDSEntity_Quad_Tetra ) continue; - if ( !volTool.Set( *elem, /*ignoreCentralNodes=*/false )) continue; // strange... + if ( !volTool.Set( elem, /*ignoreCentralNodes=*/false )) continue; // strange... - TSplitMethod splitMethod = getSplitMethod( volTool, theMethodFlags ); - if ( splitMethod._nbTetra < 1 ) continue; + TSplitMethod splitMethod = ( facetToSplit < 0 ? + getTetraSplitMethod( volTool, theMethodFlags ) : + getPrismSplitMethod( volTool, theMethodFlags, facetToSplit )); + if ( splitMethod._nbSplits < 1 ) continue; // find submesh to add new tetras to - if ( !subMesh || !subMesh->Contains( *elem )) + if ( !subMesh || !subMesh->Contains( elem )) { - int shapeID = FindShape( *elem ); + int shapeID = FindShape( elem ); helper.SetSubShape( shapeID ); // helper will add tetras to the found submesh subMesh = GetMeshDS()->MeshElements( shapeID ); } int iQ; - if ( (*elem)->IsQuadratic() ) + if ( elem->IsQuadratic() ) { iQ = 2; // add quadratic links to the helper @@ -1963,7 +2241,8 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, iQ = 1; helper.SetIsQuadratic( false ); } - vector nodes( (*elem)->begin_nodes(), (*elem)->end_nodes() ); + vector nodes( volTool.GetNodes(), + volTool.GetNodes() + elem->NbNodes() ); helper.SetElementsOnShape( true ); if ( splitMethod._baryNode ) { @@ -1991,16 +2270,25 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } } - // make tetras - vector tetras( splitMethod._nbTetra ); // splits of a volume - const int* tetConn = splitMethod._connectivity; - for ( int i = 0; i < splitMethod._nbTetra; ++i, tetConn += 4 ) - newElems.Append( tetras[ i ] = helper.AddVolume( nodes[ tetConn[0] ], - nodes[ tetConn[1] ], - nodes[ tetConn[2] ], - nodes[ tetConn[3] ])); + // make new volumes + vector splitVols( splitMethod._nbSplits ); // splits of a volume + const int* volConn = splitMethod._connectivity; + if ( splitMethod._nbCorners == 4 ) // tetra + for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners ) + newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ], + nodes[ volConn[1] ], + nodes[ volConn[2] ], + nodes[ volConn[3] ])); + else // prisms + for ( int i = 0; i < splitMethod._nbSplits; ++i, volConn += splitMethod._nbCorners ) + newElems.Append( splitVols[ i ] = helper.AddVolume( nodes[ volConn[0] ], + nodes[ volConn[1] ], + nodes[ volConn[2] ], + nodes[ volConn[3] ], + nodes[ volConn[4] ], + nodes[ volConn[5] ])); - ReplaceElemInGroups( *elem, tetras, GetMeshDS() ); + ReplaceElemInGroups( elem, splitVols, GetMeshDS() ); // Split faces on sides of the split volume @@ -2029,17 +2317,37 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, map::iterator iF_n = splitMethod._faceBaryNode.find(iF); if ( iF_n != splitMethod._faceBaryNode.end() ) { + const SMDS_MeshNode *baryNode = iF_n->second; for ( int iN = 0; iN < nbNodes*iQ; iN += iQ ) { const SMDS_MeshNode* n1 = fNodes[iN]; const SMDS_MeshNode *n2 = fNodes[(iN+iQ)%(nbNodes*iQ)]; - const SMDS_MeshNode *n3 = iF_n->second; + const SMDS_MeshNode *n3 = baryNode; if ( !volTool.IsFaceExternal( iF )) swap( n2, n3 ); triangles.push_back( helper.AddFace( n1,n2,n3 )); - - if ( fSubMesh && n3->getshapeId() < 1 ) - fSubMesh->AddNode( n3 ); + } + if ( fSubMesh ) // update position of the bary node on geometry + { + if ( subMesh ) + subMesh->RemoveNode( baryNode, false ); + GetMeshDS()->SetNodeOnFace( baryNode, fSubMesh->GetID() ); + const TopoDS_Shape& s = GetMeshDS()->IndexToShape( fSubMesh->GetID() ); + if ( !s.IsNull() && s.ShapeType() == TopAbs_FACE ) + { + fHelper.SetSubShape( s ); + gp_XY uv( 1e100, 1e100 ); + double distXYZ[4]; + if ( !fHelper.CheckNodeUV( TopoDS::Face( s ), baryNode, + uv, /*tol=*/1e-7, /*force=*/true, distXYZ ) && + uv.X() < 1e100 ) + { + // node is too far from the surface + GetMeshDS()->MoveNode( baryNode, distXYZ[1], distXYZ[2], distXYZ[3] ); + const_cast( baryNode )->SetPosition + ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() ))); + } + } } } else @@ -2069,6 +2377,8 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } } list< TTriangleFacet >::iterator facet = facets.begin(); + if ( facet == facets.end() ) + break; for ( ; facet != facets.end(); ++facet ) { if ( !volTool.IsFaceExternal( iF )) @@ -2087,11 +2397,11 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, } ReplaceElemInGroups( face, triangles, GetMeshDS() ); GetMeshDS()->RemoveFreeElement( face, fSubMesh, /*fromGroups=*/false ); - } + } // while a face based on facet nodes exists } // loop on volume faces to split them into triangles - GetMeshDS()->RemoveFreeElement( *elem, subMesh, /*fromGroups=*/false ); + GetMeshDS()->RemoveFreeElement( elem, subMesh, /*fromGroups=*/false ); if ( geomType == SMDSEntity_TriQuad_Hexa ) { @@ -2106,6 +2416,198 @@ void SMESH_MeshEditor::SplitVolumesIntoTetra (const TIDSortedElemSet & theElems, myLastCreatedElems = newElems; } +//======================================================================= +//function : GetHexaFacetsToSplit +//purpose : For hexahedra that will be split into prisms, finds facets to +// split into triangles. Only hexahedra adjacent to the one closest +// to theFacetNormal.Location() are returned. +//param [in,out] theHexas - the hexahedra +//param [in] theFacetNormal - facet normal +//param [out] theFacets - the hexahedra and found facet IDs +//======================================================================= + +void SMESH_MeshEditor::GetHexaFacetsToSplit( TIDSortedElemSet& theHexas, + const gp_Ax1& theFacetNormal, + TFacetOfElem & theFacets) +{ + #define THIS_METHOD "SMESH_MeshEditor::GetHexaFacetsToSplit(): " + + // Find a hexa closest to the location of theFacetNormal + + const SMDS_MeshElement* startHex; + { + // get SMDS_ElemIteratorPtr on theHexas + typedef const SMDS_MeshElement* TValue; + typedef TIDSortedElemSet::iterator TSetIterator; + typedef SMDS::SimpleAccessor TAccesor; + typedef SMDS_MeshElement::GeomFilter TFilter; + typedef SMDS_SetIterator < TValue, TSetIterator, TAccesor, TFilter > TElemSetIter; + SMDS_ElemIteratorPtr elemIt = SMDS_ElemIteratorPtr + ( new TElemSetIter( theHexas.begin(), + theHexas.end(), + SMDS_MeshElement::GeomFilter( SMDSGeom_HEXA ))); + + SMESH_ElementSearcher* searcher = + SMESH_MeshAlgos::GetElementSearcher( *myMesh->GetMeshDS(), elemIt ); + + startHex = searcher->FindClosestTo( theFacetNormal.Location(), SMDSAbs_Volume ); + + delete searcher; + + if ( !startHex ) + throw SALOME_Exception( THIS_METHOD "startHex not found"); + } + + // Select a facet of startHex by theFacetNormal + + SMDS_VolumeTool vTool( startHex ); + double norm[3], dot, maxDot = 0; + int facetID = -1; + for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) + if ( vTool.GetFaceNormal( iF, norm[0], norm[1], norm[2] )) + { + dot = Abs( theFacetNormal.Direction().Dot( gp_Dir( norm[0], norm[1], norm[2] ))); + if ( dot > maxDot ) + { + facetID = iF; + maxDot = dot; + } + } + if ( facetID < 0 ) + throw SALOME_Exception( THIS_METHOD "facet of startHex not found"); + + // Fill theFacets starting from facetID of startHex + + // facets used for seach of volumes adjacent to already treated ones + typedef pair< TFacetOfElem::iterator, int > TElemFacets; + typedef map< TVolumeFaceKey, TElemFacets > TFacetMap; + TFacetMap facetsToCheck; + + set facetNodes; + const SMDS_MeshElement* curHex; + + const bool allHex = ( theHexas.size() == myMesh->NbHexas() ); + + while ( startHex ) + { + // move in two directions from startHex via facetID + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + curHex = startHex; + int curFacet = facetID; + if ( is2nd ) // do not treat startHex twice + { + vTool.Set( curHex ); + if ( vTool.IsFreeFace( curFacet, &curHex )) + { + curHex = 0; + } + else + { + vTool.GetFaceNodes( curFacet, facetNodes ); + vTool.Set( curHex ); + curFacet = vTool.GetFaceIndex( facetNodes ); + } + } + while ( curHex ) + { + // store a facet to split + if ( curHex->GetGeomType() != SMDSGeom_HEXA ) + { + theFacets.insert( make_pair( curHex, -1 )); + break; + } + if ( !allHex && !theHexas.count( curHex )) + break; + + pair< TFacetOfElem::iterator, bool > facetIt2isNew = + theFacets.insert( make_pair( curHex, curFacet )); + if ( !facetIt2isNew.second ) + break; + + // remember not-to-split facets in facetsToCheck + int oppFacet = vTool.GetOppFaceIndexOfHex( curFacet ); + for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) + { + if ( iF == curFacet && iF == oppFacet ) + continue; + TVolumeFaceKey facetKey ( vTool, iF ); + TElemFacets elemFacet( facetIt2isNew.first, iF ); + pair< TFacetMap::iterator, bool > it2isnew = + facetsToCheck.insert( make_pair( facetKey, elemFacet )); + if ( !it2isnew.second ) + facetsToCheck.erase( it2isnew.first ); // adjacent hex already checked + } + // pass to a volume adjacent via oppFacet + if ( vTool.IsFreeFace( oppFacet, &curHex )) + { + curHex = 0; + } + else + { + // get a new curFacet + vTool.GetFaceNodes( oppFacet, facetNodes ); + vTool.Set( curHex ); + curFacet = vTool.GetFaceIndex( facetNodes, /*hint=*/curFacet ); + } + } + } // move in two directions from startHex via facetID + + // Find a new startHex by facetsToCheck + + startHex = 0; + facetID = -1; + TFacetMap::iterator fIt = facetsToCheck.begin(); + while ( !startHex && fIt != facetsToCheck.end() ) + { + const TElemFacets& elemFacets = fIt->second; + const SMDS_MeshElement* hex = elemFacets.first->first; + int splitFacet = elemFacets.first->second; + int lateralFacet = elemFacets.second; + facetsToCheck.erase( fIt ); + fIt = facetsToCheck.begin(); + + vTool.Set( hex ); + if ( vTool.IsFreeFace( lateralFacet, &curHex ) || + curHex->GetGeomType() != SMDSGeom_HEXA ) + continue; + if ( !allHex && !theHexas.count( curHex )) + continue; + + startHex = curHex; + + // find a facet of startHex to split + + set lateralNodes; + vTool.GetFaceNodes( lateralFacet, lateralNodes ); + vTool.GetFaceNodes( splitFacet, facetNodes ); + int oppLateralFacet = vTool.GetOppFaceIndexOfHex( lateralFacet ); + vTool.Set( startHex ); + lateralFacet = vTool.GetFaceIndex( lateralNodes, oppLateralFacet ); + + // look for a facet of startHex having common nodes with facetNodes + // but not lateralFacet + for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) + { + if ( iF == lateralFacet ) + continue; + int nbCommonNodes = 0; + const SMDS_MeshNode** nn = vTool.GetFaceNodes( iF ); + for ( int iN = 0, nbN = vTool.NbFaceNodes( iF ); iN < nbN; ++iN ) + nbCommonNodes += facetNodes.count( nn[ iN ]); + + if ( nbCommonNodes >= 2 ) + { + facetID = iF; + break; + } + } + if ( facetID < 0 ) + throw SALOME_Exception( THIS_METHOD "facet of a new startHex not found"); + } + } // while ( startHex ) +} + //======================================================================= //function : AddToSameGroups //purpose : add elemToAdd to the groups the elemInGroups belongs to @@ -3200,13 +3702,8 @@ static bool getClosestUV (Extrema_GenExtPS& projector, if ( projector.IsDone() ) { double u, v, minVal = DBL_MAX; for ( int i = projector.NbExt(); i > 0; i-- ) -#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1 if ( projector.SquareDistance( i ) < minVal ) { minVal = projector.SquareDistance( i ); -#else - if ( projector.Value( i ) < minVal ) { - minVal = projector.Value( i ); -#endif projector.Point( i ).Parameter( u, v ); } result.SetCoord( u, v ); @@ -4200,7 +4697,7 @@ void SMESH_MeshEditor::sweepElement(const SMDS_MeshElement* elem, //======================================================================= void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, - TElemOfElemListMap & newElemsMap, + TTElemOfElemListMap & newElemsMap, TElemOfVecOfNnlmiMap & elemNewNodesMap, TIDSortedElemSet& elemSet, const int nbSteps, @@ -4245,7 +4742,7 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap & mapNewNodes, // Make a ceiling for each element ie an equal element of last new nodes. // Find free links of faces - make edges and sweep them into faces. - TElemOfElemListMap::iterator itElem = newElemsMap.begin(); + TTElemOfElemListMap::iterator itElem = newElemsMap.begin(); TElemOfVecOfNnlmiMap::iterator itElemNodes = elemNewNodesMap.begin(); for ( ; itElem != newElemsMap.end(); itElem++, itElemNodes++ ) { @@ -4630,7 +5127,7 @@ SMESH_MeshEditor::RotationSweep(TIDSortedElemSet & theElems, TNodeOfNodeListMap mapNewNodes; TElemOfVecOfNnlmiMap mapElemNewNodes; - TElemOfElemListMap newElemsMap; + TTElemOfElemListMap newElemsMap; const bool isQuadraticMesh = bool( myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC) + @@ -4772,13 +5269,13 @@ const SMDS_MeshNode* SMESH_MeshEditor::CreateNode(const double x, //======================================================================= SMESH_MeshEditor::PGroupIDs -SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems, - const gp_Vec& theStep, - const int theNbSteps, - TElemOfElemListMap& newElemsMap, - const bool theMakeGroups, - const int theFlags, - const double theTolerance) +SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems, + const gp_Vec& theStep, + const int theNbSteps, + TTElemOfElemListMap& newElemsMap, + const bool theMakeGroups, + const int theFlags, + const double theTolerance) { ExtrusParam aParams; aParams.myDir = gp_Dir(theStep); @@ -4799,12 +5296,12 @@ SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems, //======================================================================= SMESH_MeshEditor::PGroupIDs -SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems, - ExtrusParam& theParams, - TElemOfElemListMap& newElemsMap, - const bool theMakeGroups, - const int theFlags, - const double theTolerance) +SMESH_MeshEditor::ExtrusionSweep (TIDSortedElemSet & theElems, + ExtrusParam& theParams, + TTElemOfElemListMap& newElemsMap, + const bool theMakeGroups, + const int theFlags, + const double theTolerance) { myLastCreatedElems.Clear(); myLastCreatedNodes.Clear(); @@ -4976,7 +5473,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, if( aS.ShapeType() == TopAbs_EDGE ) { aTrackEdge = TopoDS::Edge( aS ); // the Edge must not be degenerated - if ( BRep_Tool::Degenerated( aTrackEdge ) ) + if ( SMESH_Algo::isDegenerated( aTrackEdge ) ) return EXTR_BAD_PATH_SHAPE; TopExp::Vertices( aTrackEdge, aV1, aV2 ); aItN = theTrack->GetFather()->GetSubMesh( aV1 )->GetSubMeshDS()->GetNodes(); @@ -5264,7 +5761,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, else if( aS.ShapeType() == TopAbs_EDGE ) { aTrackEdge = TopoDS::Edge( aS ); // the Edge must not be degenerated - if ( BRep_Tool::Degenerated( aTrackEdge ) ) + if ( SMESH_Algo::isDegenerated( aTrackEdge ) ) return EXTR_BAD_PATH_SHAPE; TopExp::Vertices( aTrackEdge, aV1, aV2 ); const SMDS_MeshNode* aN1 = SMESH_Algo::VertexNode( aV1, pMeshDS ); @@ -5290,7 +5787,7 @@ SMESH_MeshEditor::ExtrusionAlongTrack (TIDSortedElemSet & theElements, TopExp_Explorer eExp(aS, TopAbs_EDGE); for(; eExp.More(); eExp.Next()) { TopoDS_Edge E = TopoDS::Edge( eExp.Current() ); - if( BRep_Tool::Degenerated(E) ) continue; + if( SMESH_Algo::isDegenerated(E) ) continue; SMESH_subMesh* SM = theTrack->GetSubMesh(E); if(SM) { LSM.push_back(SM); @@ -5463,12 +5960,12 @@ SMESH_MeshEditor::MakeExtrElements(TIDSortedElemSet& theElements for( ; itPP != fullList.end(); itPP++) { aPPs.push_back( *itPP ); if ( theHasAngles && itAngles != theAngles.end() ) - aPPs.back().SetAngle( *itAngles ); + aPPs.back().SetAngle( *itAngles++ ); } TNodeOfNodeListMap mapNewNodes; TElemOfVecOfNnlmiMap mapElemNewNodes; - TElemOfElemListMap newElemsMap; + TTElemOfElemListMap newElemsMap; TIDSortedElemSet::iterator itElem; // source elements for each generated one SMESH_SequenceOfElemPtr srcElems, srcNodes; @@ -6011,7 +6508,7 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems, if ( ( theMakeGroups && theCopy ) || ( theMakeGroups && theTargetMesh ) ) - newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh ); + newGroupIDs = generateGroups( srcNodes, srcElems, groupPostfix, theTargetMesh, false ); return newGroupIDs; } @@ -6019,9 +6516,11 @@ SMESH_MeshEditor::Transform (TIDSortedElemSet & theElems, //======================================================================= /*! * \brief Create groups of elements made during transformation - * \param nodeGens - nodes making corresponding myLastCreatedNodes - * \param elemGens - elements making corresponding myLastCreatedElems - * \param postfix - to append to names of new groups + * \param nodeGens - nodes making corresponding myLastCreatedNodes + * \param elemGens - elements making corresponding myLastCreatedElems + * \param postfix - to append to names of new groups + * \param targetMesh - mesh to create groups in + * \param topPresent - is there "top" elements that are created by sweeping */ //======================================================================= @@ -6029,14 +6528,17 @@ SMESH_MeshEditor::PGroupIDs SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, const SMESH_SequenceOfElemPtr& elemGens, const std::string& postfix, - SMESH_Mesh* targetMesh) + SMESH_Mesh* targetMesh, + const bool topPresent) { PGroupIDs newGroupIDs( new list ); SMESH_Mesh* mesh = targetMesh ? targetMesh : GetMesh(); // Sort existing groups by types and collect their names - // to store an old group and a generated new ones + // containers to store an old group and generated new ones; + // 1st new group is for result elems of different type than a source one; + // 2nd new group is for same type result elems ("top" group at extrusion) using boost::tuple; using boost::make_tuple; typedef tuple< SMESHDS_GroupBase*, SMESHDS_Group*, SMESHDS_Group* > TOldNewGroup; @@ -6066,6 +6568,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, // Loop on nodes and elements to add them in new groups + vector< const SMDS_MeshElement* > resultElems; for ( int isNodes = 0; isNodes < 2; ++isNodes ) { const SMESH_SequenceOfElemPtr& gens = isNodes ? nodeGens : elemGens; @@ -6088,7 +6591,7 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, continue; } // collect all elements made by the iElem-th sourceElem - list< const SMDS_MeshElement* > resultElems; + resultElems.clear(); if ( const SMDS_MeshElement* resElem = elems( iElem )) if ( resElem != sourceElem ) resultElems.push_back( resElem ); @@ -6097,25 +6600,23 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, if ( resElem != sourceElem ) resultElems.push_back( resElem ); - // there must be a top element const SMDS_MeshElement* topElem = 0; - if ( isNodes ) + if ( isNodes ) // there must be a top element { topElem = resultElems.back(); resultElems.pop_back(); } else { - list< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin(); + vector< const SMDS_MeshElement* >::reverse_iterator resElemIt = resultElems.rbegin(); for ( ; resElemIt != resultElems.rend() ; ++resElemIt ) if ( (*resElemIt)->GetType() == sourceElem->GetType() ) { topElem = *resElemIt; - resultElems.erase( --(resElemIt.base()) ); // erase *resElemIt + *resElemIt = 0; // erase *resElemIt break; } } - // add resultElems to groups originted from ones the sourceElem belongs to list< TOldNewGroup >::iterator gOldNew, gLast = groupsOldNew.end(); for ( gOldNew = groupsOldNew.begin(); gOldNew != gLast; ++gOldNew ) @@ -6125,16 +6626,17 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, { // fill in a new group SMDS_MeshGroup & newGroup = gOldNew->get<1>()->SMDSGroup(); - list< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt; + vector< const SMDS_MeshElement* >::iterator resLast = resultElems.end(), resElemIt; for ( resElemIt = resultElems.begin(); resElemIt != resLast; ++resElemIt ) - newGroup.Add( *resElemIt ); + if ( *resElemIt ) + newGroup.Add( *resElemIt ); // fill a "top" group if ( topElem ) { SMDS_MeshGroup & newTopGroup = gOldNew->get<2>()->SMDSGroup(); newTopGroup.Add( topElem ); - } + } } } } // loop on created elements @@ -6148,7 +6650,6 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, SMESHDS_GroupBase* oldGroupDS = orderedOldNewGroups[i]->get<0>(); SMESHDS_Group* newGroups[2] = { orderedOldNewGroups[i]->get<1>(), orderedOldNewGroups[i]->get<2>() }; - const int nbNewGroups = !newGroups[0]->IsEmpty() + !newGroups[1]->IsEmpty(); for ( int is2nd = 0; is2nd < 2; ++is2nd ) { SMESHDS_Group* newGroupDS = newGroups[ is2nd ]; @@ -6162,11 +6663,21 @@ SMESH_MeshEditor::generateGroups(const SMESH_SequenceOfElemPtr& nodeGens, newGroupDS->SetType( newGroupDS->GetElements()->next()->GetType() ); // make a name - const bool isTop = ( nbNewGroups == 2 && + const bool isTop = ( topPresent && newGroupDS->GetType() == oldGroupDS->GetType() && is2nd ); string name = oldGroupDS->GetStoreName(); + { // remove trailing whitespaces (issue 22599) + size_t size = name.size(); + while ( size > 1 && isspace( name[ size-1 ])) + --size; + if ( size != name.size() ) + { + name.resize( size ); + oldGroupDS->SetStoreName( name.c_str() ); + } + } if ( !targetMesh ) { string suffix = ( isTop ? "top": postfix.c_str() ); name += "_"; @@ -8334,6 +8845,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT aHelper.SetIsQuadratic( true ); aHelper.SetIsBiQuadratic( theToBiQuad ); aHelper.SetElementsOnShape(true); + aHelper.ToFixNodeParameters( true ); // convert elements assigned to sub-meshes int nbCheckedElems = 0; @@ -8446,11 +8958,20 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theT if ( !volume ) continue; const SMDSAbs_EntityType type = volume->GetEntityType(); - if (( theToBiQuad && type == SMDSEntity_TriQuad_Hexa ) || - ( !theToBiQuad && type == SMDSEntity_Quad_Hexa )) + if ( volume->IsQuadratic() ) { - aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume )); - continue; + bool alreadyOK; + switch ( type ) + { + case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break; + case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break; + default: alreadyOK = true; + } + if ( alreadyOK ) + { + aHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( volume )); + continue; + } } const int id = volume->GetID(); vector nodes (volume->begin_nodes(), volume->end_nodes()); @@ -9708,6 +10229,7 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh* theMeshDS, { // duplicate node aNewNode = theMeshDS->AddNode( aCurrNode->X(), aCurrNode->Y(), aCurrNode->Z() ); + copyPosition( aCurrNode, aNewNode ); theNodeNodeMap[ aCurrNode ] = aNewNode; myLastCreatedNodes.Append( aNewNode ); } @@ -9720,10 +10242,8 @@ bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh* theMeshDS, if ( theIsDoubleElem ) AddElement(newNodes, anElem->GetType(), anElem->IsPoly()); else - { - MESSAGE("ChangeElementNodes"); theMeshDS->ChangeElementNodes( anElem, &newNodes[ 0 ], anElem->NbNodes() ); - } + res = true; } return res; @@ -9734,8 +10254,8 @@ 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 - they not assigned to elements + 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 */ //================================================================================ @@ -9771,6 +10291,7 @@ bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes, const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() ); if ( aNewNode ) { + copyPosition( aNode, aNewNode ); anOldNodeToNewNode[ aNode ] = aNewNode; myLastCreatedNodes.Append( aNewNode ); } @@ -9872,15 +10393,12 @@ namespace { } void Perform(const gp_Pnt& aPnt, double theTol) { + theTol *= theTol; _state = TopAbs_OUT; _extremum.Perform(aPnt); if ( _extremum.IsDone() ) for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol) -#if OCC_VERSION_LARGE > 0x06040000 // Porting to OCCT6.5.1 _state = ( _extremum.SquareDistance(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT ); -#else - _state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT ); -#endif } TopAbs_State State() const { @@ -9891,12 +10409,14 @@ namespace { //================================================================================ /*! - \brief Identify the elements that will be affected by node duplication (actual duplication is not performed. + \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). + located on or inside shape). If the shape is null, detection is done on faces orientations + (select elements with a gravity center on the side given by faces normals). + This mode (null shape) is faster, but works only when theElems are faces, with coherents orientations. The replicated nodes should be associated to affected elements. \return groups of affected elements \sa DoubleNodeElemGroupsInRegion() @@ -9909,44 +10429,145 @@ bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theEl 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))); + std::set alreadyCheckedNodes; + std::set alreadyCheckedElems; + std::set edgesToCheck; + alreadyCheckedNodes.clear(); + alreadyCheckedElems.clear(); + edgesToCheck.clear(); + + // --- iterates on elements to be replicated and get elements by back references from their nodes + + TIDSortedElemSet::const_iterator elemItr = theElems.begin(); + int ielem; + for ( ielem=1; elemItr != theElems.end(); ++elemItr ) + { + SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr; + if (!anElem || (anElem->GetType() != SMDSAbs_Face)) + continue; + gp_XYZ normal; + SMESH_MeshAlgos::FaceNormal( anElem, normal, /*normalized=*/true ); + MESSAGE("element " << ielem++ << " normal " << normal.X() << " " << normal.Y() << " " << normal.Z()); + std::set nodesElem; + nodesElem.clear(); + SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator(); + while ( nodeItr->more() ) + { + const SMDS_MeshNode* aNode = cast2Node(nodeItr->next()); + nodesElem.insert(aNode); + } + std::set::iterator nodit = nodesElem.begin(); + for (; nodit != nodesElem.end(); nodit++) + { + MESSAGE(" noeud "); + const SMDS_MeshNode* aNode = *nodit; + if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() ) + continue; + if (alreadyCheckedNodes.find(aNode) != alreadyCheckedNodes.end()) + continue; + alreadyCheckedNodes.insert(aNode); + SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator(); + while ( backElemItr->more() ) + { + MESSAGE(" backelem "); + const SMDS_MeshElement* curElem = backElemItr->next(); + if (alreadyCheckedElems.find(curElem) != alreadyCheckedElems.end()) + continue; + if (theElems.find(curElem) != theElems.end()) + continue; + alreadyCheckedElems.insert(curElem); + double x=0, y=0, z=0; + int nb = 0; + SMDS_ElemIteratorPtr nodeItr2 = curElem->nodesIterator(); + while ( nodeItr2->more() ) + { + const SMDS_MeshNode* anotherNode = cast2Node(nodeItr2->next()); + x += anotherNode->X(); + y += anotherNode->Y(); + z += anotherNode->Z(); + nb++; + } + gp_XYZ p; + p.SetCoord( x/nb -aNode->X(), + y/nb -aNode->Y(), + z/nb -aNode->Z() ); + MESSAGE(" check " << p.X() << " " << p.Y() << " " << p.Z()); + if (normal*p > 0) + { + MESSAGE(" --- inserted") + theAffectedElems.insert( curElem ); + } + else if (curElem->GetType() == SMDSAbs_Edge) + edgesToCheck.insert(curElem); + } + } + } + // --- add also edges lying on the set of faces (all nodes in alreadyCheckedNodes) + std::set::iterator eit = edgesToCheck.begin(); + for( ; eit != edgesToCheck.end(); eit++) + { + bool onside = true; + const SMDS_MeshElement* anEdge = *eit; + SMDS_ElemIteratorPtr nodeItr = anEdge->nodesIterator(); + while ( nodeItr->more() ) + { + const SMDS_MeshNode* aNode = cast2Node(nodeItr->next()); + if (alreadyCheckedNodes.find(aNode) == alreadyCheckedNodes.end()) + { + onside = false; + break; + } + } + if (onside) + { + MESSAGE(" --- edge onside inserted") + theAffectedElems.insert(anEdge); + } + } } - - // iterates on indicated elements and get elements by back references from their nodes - TIDSortedElemSet::const_iterator elemItr = theElems.begin(); - for ( ; elemItr != theElems.end(); ++elemItr ) + else { - SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr; - if (!anElem) - continue; + 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))); + } - SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator(); - while ( nodeItr->more() ) + // iterates on indicated elements and get elements by back references from their nodes + TIDSortedElemSet::const_iterator elemItr = theElems.begin(); + int ielem; + for ( ielem = 1; elemItr != theElems.end(); ++elemItr ) { - const SMDS_MeshNode* aNode = cast2Node(nodeItr->next()); - if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() ) + MESSAGE("element " << ielem++); + SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr; + if (!anElem) continue; - SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator(); - while ( backElemItr->more() ) + SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator(); + while ( nodeItr->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 ); + MESSAGE(" noeud "); + const SMDS_MeshNode* aNode = cast2Node(nodeItr->next()); + if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() ) + continue; + SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator(); + while ( backElemItr->more() ) + { + MESSAGE(" backelem "); + 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 ); + } } } } @@ -10034,7 +10655,12 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const gp_Vec v2(p0, g2); gp_Vec n1 = vref.Crossed(v1); gp_Vec n2 = vref.Crossed(v2); - return n2.AngleWithRef(n1, vref); + try { + return n2.AngleWithRef(n1, vref); + } + catch ( Standard_Failure ) { + } + return Max( v1.Magnitude(), v2.Magnitude() ); } /*! @@ -10047,13 +10673,16 @@ double SMESH_MeshEditor::OrientedAngle(const gp_Pnt& p0, const gp_Pnt& p1, const * If there is no shared faces between the group #n and the group #p in the list, the group j_n_p is not created. * All the flat elements are gathered into the group named "joints3D" (or "joints2D" in 2D situation). * The flat element of the multiple junctions between the simple junction are stored in a group named "jointsMultiples". - * @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 - * @return TRUE if operation has been completed successfully, FALSE otherwise + * \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 + * \param onAllBoundaries - if TRUE, the nodes and elements are also created on + * the boundary between \a theDomains and the rest mesh + * \return TRUE if operation has been completed successfully, FALSE otherwise */ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector& theElems, - bool createJointElems) + bool createJointElems, + bool onAllBoundaries) { MESSAGE("----------------------------------------------"); MESSAGE("SMESH_MeshEditor::doubleNodesOnGroupBoundaries"); @@ -10082,15 +10711,20 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector volume to modify) @@ -10119,7 +10753,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorgetVtkId(); @@ -10132,26 +10766,30 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorfromVtkToSmds(neighborsVtkIds[n]); const SMDS_MeshElement* elem = meshDS->FindElement(smdsId); - if (! domain.count(elem)) // neighbor is in another domain : face is shared + if (elem && ! domain.count(elem)) // neighbor is in another domain : face is shared { bool ok = false ; - for (int idombis = 0; idombis < theElems.size(); idombis++) // check if the neighbor belongs to another domain of the list + for (int idombis = 0; idombis < theElems.size() && !ok; idombis++) // check if the neighbor belongs to another domain of the list { // MESSAGE("Domain " << idombis); const TIDSortedElemSet& domainbis = theElems[idombis]; if ( domainbis.count(elem)) ok = true ; // neighbor is in a correct domain : face is kept } - if ( ok ) // the characteristics of the face is stored + if ( ok || onAllBoundaries ) // the characteristics of the face is stored { DownIdType face(downIds[n], downTypes[n]); - if (!faceDomains.count(face)) - faceDomains[face] = emptyMap; // create an empty entry for face if (!faceDomains[face].count(idom)) { faceDomains[face][idom] = vtkId; // volume associated to face in this domain celldom[vtkId] = idom; //MESSAGE(" cell with a border " << vtkId << " domain " << idom); } + if ( !ok ) + { + theRestDomElems.insert( elem ); + faceDomains[face][iRestDom] = neighborsVtkIds[n]; + celldom[neighborsVtkIds[n]] = iRestDom; + } } } } @@ -10165,14 +10803,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector domvol = itface->second; + const std::map& domvol = itface->second; if (!domvol.count(idomain)) continue; DownIdType face = itface->first; @@ -10201,8 +10839,6 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector > mutipleNodesToFace; // nodes multi domains with domain order to transform in Face (junction between 3 or more 2D domains) MESSAGE(".. Duplication of the nodes"); - for (int idomain = 0; idomain < theElems.size(); idomain++) + for (int idomain = idom0; idomain < nbDomains; idomain++) { itface = faceDomains.begin(); for (; itface != faceDomains.end(); ++itface) { - std::map domvol = itface->second; + const std::map& domvol = itface->second; if (!domvol.count(idomain)) continue; DownIdType face = itface->first; @@ -10241,15 +10877,12 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector::iterator itdom = domvol.begin(); + std::map::const_iterator itdom = domvol.begin(); for (; itdom != domvol.end(); ++itdom) { int idom = itdom->first; @@ -10277,6 +10910,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vectorGetPoint(oldId); SMDS_MeshNode *newNode = meshDS->AddNode(coords[0], coords[1], coords[2]); + copyPosition( meshDS->FindNodeVtk( oldId ), newNode ); int newId = newNode->getVtkId(); nodeDomains[oldId][idom] = newId; // cloned node for other domains //MESSAGE("-+-+-c oldNode " << oldId << " domain " << idomain << " newNode " << newId << " domain " << idom << " size=" <fromVtkToSmds(vtkVolIds[ivol]); SMDS_MeshElement* elem = (SMDS_MeshElement*)meshDS->FindElement(smdsId); - if (theElems[idom].count(elem)) + if (domain.count(elem)) { SMDS_VtkVolume* svol = dynamic_cast(elem); domvol[idom] = svol; @@ -10541,7 +11176,7 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector >::const_iterator itnod = nodeDomains.begin(); for (; itnod != nodeDomains.end(); ++itnod) @@ -10614,6 +11249,14 @@ bool SMESH_MeshEditor::DoubleNodesOnGroupBoundaries( const std::vector::iterator name_group = mapOfJunctionGroups.begin(); + for ( ; name_group != mapOfJunctionGroups.end(); ++name_group ) + { + if ( name_group->second && name_group->second->GetGroupDS()->IsEmpty() ) + myMesh->RemoveGroup( name_group->second->GetGroupDS()->GetID() ); + } + meshDS->CleanDownWardConnectivity(); // Mesh has been modified, downward connectivity is no more usable, free memory grid->BuildLinks(); @@ -10681,6 +11324,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vectorAddNode(node->X(), node->Y(), node->Z()); + copyPosition( node, clone ); clonedNodes[node] = clone; } else @@ -10697,6 +11341,7 @@ bool SMESH_MeshEditor::CreateFlatElementsOnFacesGroups(const std::vectorAddNode(node->X(), node->Y(), node->Z()); + copyPosition( node, inter ); intermediateNodes[node] = inter; } else @@ -11641,3 +12286,45 @@ int SMESH_MeshEditor::MakeBoundaryMesh(const TIDSortedElemSet& elements, } return nbAddedBnd; } + +//================================================================================ +/*! + * \brief Copy node position and set \a to node on the same geometry + */ +//================================================================================ + +void SMESH_MeshEditor::copyPosition( const SMDS_MeshNode* from, + const SMDS_MeshNode* to ) +{ + if ( !from || !to ) return; + + SMDS_PositionPtr pos = from->GetPosition(); + if ( !pos || from->getshapeId() < 1 ) return; + + switch ( pos->GetTypeOfPosition() ) + { + case SMDS_TOP_3DSPACE: break; + + case SMDS_TOP_FACE: + { + const SMDS_FacePosition* fPos = static_cast< const SMDS_FacePosition* >( pos ); + GetMeshDS()->SetNodeOnFace( to, from->getshapeId(), + fPos->GetUParameter(), fPos->GetVParameter() ); + break; + } + case SMDS_TOP_EDGE: + { + // WARNING: it is dangerous to set equal nodes on one EDGE!!!!!!!! + const SMDS_EdgePosition* ePos = static_cast< const SMDS_EdgePosition* >( pos ); + GetMeshDS()->SetNodeOnEdge( to, from->getshapeId(), ePos->GetUParameter() ); + break; + } + case SMDS_TOP_VERTEX: + { + GetMeshDS()->SetNodeOnVertex( to, from->getshapeId() ); + break; + } + case SMDS_TOP_UNSPEC: + default:; + } +}