X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_HexaFromSkin_3D.cxx;h=6e8448e0a2092d0bf0226dc03dc77a8393f55df2;hp=28b628c5c3d476f98272092ea12d5be06effe1a0;hb=0c16e57723ec8522693f2f0b6d3216b4f4e3d686;hpb=a399364c85efef719b8f0f1d1a81fc57a60b4601 diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx index 28b628c5c..6e8448e0a 100644 --- a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx +++ b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx @@ -51,7 +51,7 @@ namespace { enum EBoxSides //!< sides of the block { - B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED + B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES }; const char* SBoxSides[] = //!< names of block sides { @@ -59,7 +59,7 @@ namespace }; enum EQuadEdge //!< edges of quadrangle side { - Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED + Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES }; @@ -84,9 +84,42 @@ namespace return true; } + //================================================================================ + /*! + * \brief Description of node used to detect corner nodes + */ + struct _NodeDescriptor + { + int nbInverseFaces, nbNodesInInverseFaces; + + _NodeDescriptor(const SMDS_MeshNode* n): nbInverseFaces(0), nbNodesInInverseFaces(0) + { + if ( n ) + { + set nodes; + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + nodes.insert( face->begin_nodes(), face->end_nodes() ); + nbInverseFaces++; + } + nbNodesInInverseFaces = nodes.size(); + } + } + bool operator==(const _NodeDescriptor& other) const + { + return + nbInverseFaces == other.nbInverseFaces && + nbNodesInInverseFaces == other.nbNodesInInverseFaces; + } + }; + //================================================================================ /*! * \brief return true if a node is at block corner + * + * This check is valid for simple cases only */ //================================================================================ @@ -103,7 +136,7 @@ namespace bool isQuadrangle(const SMDS_MeshElement* e) { - return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 )); + return ( e && e->NbCornerNodes() == 4 ); } //================================================================================ @@ -189,15 +222,15 @@ namespace int _nbBlocksExpected; int _nbBlocksFound; -#ifdef _DEBUG_ -#define _grid_access_(args) _grid.at( args ) +#ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index +#define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)] #else -#define _grid_access_(args) _grid[ args ] +#define _grid_access_(pobj, i) pobj->_grid[ i ] #endif //!< Return node at XY - const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); } + const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));} //!< Set node at XY - void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; } + void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; } //!< Return an edge SMESH_OrientedLink getEdge(EQuadEdge edge) const { @@ -242,7 +275,7 @@ namespace //!< return coordinates by XY gp_XYZ xyz(int x, int y) const { - return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) ); + return SMESH_MeshEditor::TNodeXYZ( _grid_access_(_side, _index( x, y )) ); } //!< safely return a node by XY const SMDS_MeshNode* node(int x, int y) const @@ -259,7 +292,7 @@ namespace //!< Return a corner node const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const { - return _side->_grid_access_( _index.corner( isXMax, isYMax )); + return _grid_access_(_side, _index.corner( isXMax, isYMax )); } //!< return its size in nodes int getHoriSize() const { return _index.xSize(); } @@ -276,14 +309,29 @@ namespace */ struct _Block { - _OrientedBlockSide _side[6]; // 6 sides of a sub-block + _OrientedBlockSide _side[6]; // 6 sides of a sub-block + set _corners; const _OrientedBlockSide& getSide(int i) const { return _side[i]; } + bool setSide( int i, const _OrientedBlockSide& s) + { + if (( _side[i] = s )) + { + _corners.insert( s.cornerNode(0,0)); + _corners.insert( s.cornerNode(1,0)); + _corners.insert( s.cornerNode(0,1)); + _corners.insert( s.cornerNode(1,1)); + } + return s; + } + void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); } bool hasSide( const _OrientedBlockSide& s) const { if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true; return false; } + int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; } + bool isValid() const; }; //================================================================================ /*! @@ -300,7 +348,9 @@ namespace const SMESH_Comment& error() const { return _error; } private: - bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad); + bool fillSide( _BlockSide& side, + const SMDS_MeshElement* cornerQuad, + const SMDS_MeshNode* cornerNode); bool fillRowsUntilCorner(const SMDS_MeshElement* quad, const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, @@ -309,16 +359,14 @@ namespace bool alongN1N2 ); _OrientedBlockSide findBlockSide( EBoxSides startBlockSide, EQuadEdge sharedSideEdge1, - EQuadEdge sharedSideEdge2); + EQuadEdge sharedSideEdge2, + bool withGeometricAnalysis); //!< update own data and data of the side bound to block void setSideBoundToBlock( _BlockSide& side ) { - side._nbBlocksFound++; - if ( side.isBound() ) - { - for ( int e = 0; e < int(Q_UNDEFINED); ++e ) + if ( side._nbBlocksFound++, side.isBound() ) + for ( int e = 0; e < int(NB_QUAD_SIDES); ++e ) _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side ); - } } //!< store reason of error int error(const SMESH_Comment& reason) { _error = reason; return 0; } @@ -383,7 +431,7 @@ namespace _allSides.push_back( _BlockSide() ); _BlockSide& side = _allSides.back(); - if ( !fillSide( side, face ) ) + if ( !fillSide( side, face, *corner ) ) { if ( !_error.empty() ) return false; @@ -394,13 +442,10 @@ namespace for ( int isYMax = 0; isYMax < 2; ++isYMax ) { const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax ); - if ( !isCornerNode( nCorner )) - return BAD_MESH_ERR; - //_corner2sides[ nCorner ].insert( &side ); corners.push_back( nCorner ); cornerFaces.insert( side.getCornerFace( nCorner )); } - for ( int e = 0; e < int(Q_UNDEFINED); ++e ) + for ( int e = 0; e < int(NB_QUAD_SIDES); ++e ) _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side ); nbFacesOnSides += side.getNbFaces(); @@ -425,23 +470,33 @@ namespace return BAD_MESH_ERR; if ( _allSides.back()._grid.empty() ) _allSides.pop_back(); + _DUMP_("Nb detected sides "<< _allSides.size()); // --------------------------- // Organize sides into blocks // --------------------------- - // analyse sharing of sides by blocks - int nbBlockSides = 0; // nb of block sides taking into account their sharing - list < _BlockSide >::iterator sideIt = _allSides.begin(); - for ( ; sideIt != _allSides.end(); ++sideIt ) + // analyse sharing of sides by blocks and sort sides by nb of adjacent sides + int nbBlockSides = 0; // total nb of block sides taking into account their sharing + multimap sortedSides; { - _BlockSide& side = *sideIt; - bool isSharedSide = true; - for ( int e = 0; e < int(Q_UNDEFINED) && isSharedSide; ++e ) - isSharedSide = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size() > 2; - side._nbBlocksFound = 0; - side._nbBlocksExpected = isSharedSide ? 2 : 1; - nbBlockSides += side._nbBlocksExpected; + list < _BlockSide >::iterator sideIt = _allSides.begin(); + for ( ; sideIt != _allSides.end(); ++sideIt ) + { + _BlockSide& side = *sideIt; + bool isSharedSide = true; + int nbAdjacent = 0; + for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e ) + { + int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size(); + nbAdjacent += nbAdj; + isSharedSide = ( nbAdj > 2 ); + } + side._nbBlocksFound = 0; + side._nbBlocksExpected = isSharedSide ? 2 : 1; + nbBlockSides += side._nbBlocksExpected; + sortedSides.insert( make_pair( nbAdjacent, & side )); + } } // find sides of each block @@ -449,9 +504,9 @@ namespace while ( nbBlockSides >= 6 ) { // get any side not bound to all blocks it belongs to - sideIt = _allSides.begin(); - while ( sideIt != _allSides.end() && sideIt->isBound()) - ++sideIt; + multimap::iterator i_side = sortedSides.begin(); + while ( i_side != sortedSides.end() && i_side->second->isBound()) + ++i_side; // start searching for block sides from the got side bool ok = true; @@ -459,37 +514,40 @@ namespace _blocks.resize( _blocks.size() + 1 ); _Block& block = _blocks.back(); - block._side[B_FRONT] = &(*sideIt); - setSideBoundToBlock( *sideIt ); + block.setSide( B_FRONT, i_side->second ); + setSideBoundToBlock( *i_side->second ); nbBlockSides--; - - // edges of neighbour sides of B_FRONT corresponding to front's edges - EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT }; - EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT }; - for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i ) - ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i])); - if ( ok ) - ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP )); + // edges of adjacent sides of B_FRONT corresponding to front's edges + EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT }; + EQuadEdge edgeOfAdj [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT }; + // first find all sides detectable w/o advanced analysis, + // then repeat the search, which then may pass without advanced analysis + for ( int advAnalys = 0; advAnalys < 2; ++advAnalys ) + { + for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i ) + if ( !block._side[i] ) // try to find 4 sides adjacent to front side + ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i],edgeOfAdj[i],advAnalys)); + if ( ok || !advAnalys) + if ( !block._side[B_BACK] && block._side[B_TOP] ) // try to find back side by top one + ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP, advAnalys )); + if ( !advAnalys ) ok = true; + } + ok = block.isValid(); if ( ok ) { // check if just found block is same as one of previously found bool isSame = false; for ( int i = 1; i < _blocks.size() && !isSame; ++i ) - { - _Block& prevBlock = _blocks[i-1]; - isSame = true; - for ( int j = 0; j < 6 && isSame; ++j ) - isSame = prevBlock.hasSide( block._side[ j ]); - } + isSame = ( block._corners == _blocks[i-1]._corners ); ok = !isSame; } // count the found sides - _DUMP_(endl); - for (int i = 0; i < B_UNDEFINED; ++i ) + _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid()); + for (int i = 0; i < NB_BLOCK_SIDES; ++i ) { - _DUMP_("** Block side "<< SBoxSides[i] <<" "<< block._side[ i ]._side); + _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side); if ( block._side[ i ] ) { if ( ok && i != B_FRONT) @@ -497,22 +555,24 @@ namespace setSideBoundToBlock( *block._side[ i ]._side ); nbBlockSides--; } - _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0)); - _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0)); - _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1)); - _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1)); + _DUMP_("\t corners "<< + block._side[ i ].cornerNode(0,0)->GetID() << ", " << + block._side[ i ].cornerNode(1,0)->GetID() << ", " << + block._side[ i ].cornerNode(1,1)->GetID() << ", " << + block._side[ i ].cornerNode(0,1)->GetID() << ", "<GetNode(i))) - nCorner = firstQuad->GetNode(i); - if ( !nCorner ) return false; - // get a node on block edge int iCorner = firstQuad->GetNodeIndex( nCorner ); const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4); @@ -611,7 +666,14 @@ namespace } } - return true; + // check side validity + bool ok = + side.getCornerFace( side.getCornerNode( 0, 0 )) && + side.getCornerFace( side.getCornerNode( 1, 0 )) && + side.getCornerFace( side.getCornerNode( 0, 1 )) && + side.getCornerFace( side.getCornerNode( 1, 1 )); + + return ok; } //================================================================================ @@ -622,7 +684,8 @@ namespace _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide, EQuadEdge sharedSideEdge1, - EQuadEdge sharedSideEdge2) + EQuadEdge sharedSideEdge2, + bool withGeometricAnalysis) { _Block& block = _blocks.back(); _OrientedBlockSide& side1 = block._side[ startBlockSide ]; @@ -637,67 +700,93 @@ namespace set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set // exclude loaded sides of block from sidesOnEdge - int nbLoadedSides = 0; - for (int i = 0; i < B_UNDEFINED; ++i ) - { + for (int i = 0; i < NB_BLOCK_SIDES; ++i ) if ( block._side[ i ] ) - { - nbLoadedSides++; sidesOnEdge.erase( block._side[ i ]._side ); - } - } + int nbSidesOnEdge = sidesOnEdge.size(); _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() ); if ( nbSidesOnEdge == 0 ) return 0; _BlockSide* foundSide = 0; - if ( nbSidesOnEdge == 1 /*|| nbSidesOnEdge == 2 && nbLoadedSides == 1 */) + if ( nbSidesOnEdge == 1 ) { foundSide = *sidesOnEdge.begin(); } else { - // Select one of found sides most close to startBlockSide - - // gravity center of already loaded block sides - gp_XYZ gc(0,0,0); - for (int i = 0; i < B_UNDEFINED; ++i ) - if ( block._side[ i ] ) - gc += block._side[ i ]._side->getGC(); - gc /= nbLoadedSides; - - gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z()); - gp_Vec p1p2( p1, p2 ); - - const SMDS_MeshElement* face1 = side1->getCornerFace( n1 ); - gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1))); - gp_Vec side1Dir( p1, p1Op ); - gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir - _DUMP_(" Select adjacent for "<< side1._side << " - side dir (" - << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" ); - - map < double , _BlockSide* > angleOfSide; set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin(); - for (; sideIt != sidesOnEdge.end(); ++sideIt ) + int nbLoadedSides = block.nbSides(); + if ( nbLoadedSides > 1 ) { - _BlockSide* sideI = *sideIt; - const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 ); - gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1))); - gp_Vec sideIDir( p1, p1Op ); - // compute angle of (sideIDir projection to pln) and (X dir of pln) - gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection()); - double angle = sideIDirProj.Angle( gp::DX2d() ); - if ( angle < 0 ) angle += 2 * PI; - angleOfSide.insert( make_pair( angle, sideI )); - _DUMP_(" "<< sideI << " - side dir (" - << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")" - << " angle " << angle); + // Find the side having more than 2 corners common with already loaded sides + for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt ) + { + _BlockSide* sideI = *sideIt; + int nbCommonCorners = + block._corners.count( sideI->getCornerNode(0,0)) + + block._corners.count( sideI->getCornerNode(1,0)) + + block._corners.count( sideI->getCornerNode(0,1)) + + block._corners.count( sideI->getCornerNode(1,1)); + if ( nbCommonCorners > 2 ) + foundSide = sideI; + } + } + + if ( !foundSide ) + { + if ( !withGeometricAnalysis ) return 0; + + // Select one of found sides most close to startBlockSide + + gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z()); + gp_Vec p1p2( p1, p2 ); + + const SMDS_MeshElement* face1 = side1->getCornerFace( n1 ); + gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1))); + gp_Vec side1Dir( p1, p1Op ); + gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir + _DUMP_(" Select adjacent for "<< side1._side << " - side dir (" + << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" ); + + map < double , _BlockSide* > angleOfSide; + for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt ) + { + _BlockSide* sideI = *sideIt; + const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 ); + gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1))); + gp_Vec sideIDir( p1, p1Op ); + // compute angle of (sideIDir projection to pln) and (X dir of pln) + gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection()); + double angle = sideIDirProj.Angle( gp::DX2d() ); + if ( angle < 0 ) angle += 2 * PI; // angle [0-2*PI] + angleOfSide.insert( make_pair( angle, sideI )); + _DUMP_(" "<< sideI << " - side dir (" + << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")" + << " angle " << angle); + } + if ( nbLoadedSides == 1 ) + { + double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first; + if ( angF > PI ) angF = 2*PI - angF; + if ( angL > PI ) angL = 2*PI - angL; + foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second; + } + else + { + gp_XYZ gc(0,0,0); // gravity center of already loaded block sides + for (int i = 0; i < NB_BLOCK_SIDES; ++i ) + if ( block._side[ i ] ) + gc += block._side[ i ]._side->getGC(); + gc /= nbLoadedSides; + + gp_Vec gcDir( p1, gc ); + gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection()); + double gcAngle = gcDirProj.Angle( gp::DX2d() ); + foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second; + } } - gp_Vec gcDir( p1, gc ); - gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection()); - double gcAngle = gcDirProj.Angle( gp::DX2d() ); - foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second; _DUMP_(" selected "<< foundSide ); } @@ -759,9 +848,14 @@ namespace row2.push_back( n1 = oppositeNode( quad, i1 )); } + _NodeDescriptor nDesc( row1[1] ); + if ( nDesc == _NodeDescriptor( row1[0] )) + return true; // row of 2 nodes + // Find the rest nodes TIDSortedElemSet emptySet, avoidSet; - while ( !isCornerNode( n2 )) + //while ( !isCornerNode( n2 )) + while ( nDesc == _NodeDescriptor( n2 )) { avoidSet.clear(); avoidSet.insert( quad ); quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 ); @@ -804,7 +898,31 @@ namespace return SMDS_Mesh::FindFace(n1, n2, n3, n4 ); } -} + //================================================================================ + /*! + * \brief Checks own validity + */ + //================================================================================ + + bool _Block::isValid() const + { + bool ok = ( nbSides() == 6 ); + + // check only corners depending on side selection + EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT }; + EQuadEdge edgeAdj [4] = { Q_TOP, Q_RIGHT, Q_TOP, Q_RIGHT }; + EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT }; + + for ( int i=0; ok && i < NB_QUAD_SIDES; ++i ) + { + SMESH_OrientedLink eBack = _side[ B_BACK ].edge( edgeBack[i] ); + SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] ); + ok = ( eBack == eAdja ); + } + return ok; + } + +} // namespace //======================================================================= //function : StdMeshers_HexaFromSkin_3D