X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_HexaFromSkin_3D.cxx;h=58a90931c73d32924e057b65e8d2cf16bf66de48;hp=f0bd513ca2996ace2290e8f22cd39f5ed1a3a12d;hb=d4a710ce52f6e76786a7b3845e2f7975dc9a00b1;hpb=b87176dacc750a61d25b21f3f5b7b0d8391e07b1 diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx index f0bd513ca..58a90931c 100644 --- a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx +++ b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx @@ -1,26 +1,26 @@ -// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE // -// 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. +// 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. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // // File : StdMeshers_HexaFromSkin_3D.cxx // Created : Wed Jan 27 12:28:07 2010 // Author : Edward AGAPOV (eap) -// + #include "StdMeshers_HexaFromSkin_3D.hxx" #include "SMDS_VolumeOfNodes.hxx" @@ -33,19 +33,24 @@ //#include "utilities.h" #include -// Define error message +// Define error message and _MYDEBUG_ if needed #ifdef _DEBUG_ #define BAD_MESH_ERR \ error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \ __FILE__ ":" )<<__LINE__) +//#define _MYDEBUG_ #else #define BAD_MESH_ERR \ error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh")) #endif -// Debug output -#define _DUMP_(msg) //cout << msg << endl +// Debug output +#ifdef _MYDEBUG_ +#define _DUMP_(msg) cout << msg << endl +#else +#define _DUMP_(msg) +#endif namespace @@ -54,10 +59,12 @@ namespace { B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES }; -// const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only -// { -// "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED" -// }; +#ifdef _MYDEBUG_ + const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only + { + "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED" + }; +#endif enum EQuadEdge //!< edges of quadrangle side { Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES @@ -85,37 +92,6 @@ 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 @@ -126,7 +102,19 @@ namespace bool isCornerNode( const SMDS_MeshNode* n ) { - return n && n->NbInverseElements( SMDSAbs_Face ) % 2; + int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1; + if ( nbF % 2 ) + return true; + + set nodesInInverseFaces; + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face ); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() ); + } + + return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 ); } //================================================================================ @@ -208,9 +196,9 @@ namespace typedef void (*TFun)(int& x, int& y); TFun _xRevFun, _yRevFun, _swapFun; - static void lazy(int&, int&) {} + static void lazy (int&, int&) {} static void reverse(int& x, int& size) { x = size - x - 1; } - static void swap(int& x, int& y) { std::swap(x,y); } + static void swap (int& x, int& y) { std::swap(x,y); } }; //================================================================================ /*! @@ -244,10 +232,10 @@ namespace return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 ); } const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const; - //!< True if all blocks this side belongs to have beem found + //!< True if all blocks this side belongs to have been found bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; } //!< Return coordinates of node at XY - gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); } + gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); } //!< Return gravity center of the four corners and the middle node gp_XYZ getGC() const { @@ -276,7 +264,7 @@ namespace //!< return coordinates by XY gp_XYZ xyz(int x, int y) const { - return SMESH_MeshEditor::TNodeXYZ( _grid_access_(_side, _index( x, y )) ); + return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) ); } //!< safely return a node by XY const SMDS_MeshNode* node(int x, int y) const @@ -358,10 +346,11 @@ namespace vector& verRow1, vector& verRow2, bool alongN1N2 ); - _OrientedBlockSide findBlockSide( EBoxSides startBlockSide, - EQuadEdge sharedSideEdge1, - EQuadEdge sharedSideEdge2, - bool withGeometricAnalysis); + _OrientedBlockSide findBlockSide( EBoxSides startBlockSide, + EQuadEdge sharedSideEdge1, + EQuadEdge sharedSideEdge2, + bool withGeometricAnalysis, + set< _BlockSide* >& sidesAround); //!< update own data and data of the side bound to block void setSideBoundToBlock( _BlockSide& side ) { @@ -524,20 +513,25 @@ namespace 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 + set< _BlockSide* > sidesAround; for ( int advAnalys = 0; advAnalys < 2; ++advAnalys ) { + // try to find 4 sides adjacent to a FRONT side 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 ( !block._side[i] ) + ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i], + advAnalys, sidesAround)); + // try to find a BACK side by a TOP one 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 ( !block._side[B_BACK] && block._side[B_TOP] ) + ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP, + advAnalys, sidesAround )); if ( !advAnalys ) ok = true; } ok = block.isValid(); if ( ok ) { - // check if just found block is same as one of previously found + // check if just found block is same as one of previously found blocks bool isSame = false; for ( int i = 1; i < _blocks.size() && !isSame; ++i ) isSame = ( block._corners == _blocks[i-1]._corners ); @@ -677,16 +671,77 @@ namespace return ok; } + //================================================================================ + /*! + * \brief Return true if it's possible to make a loop over corner2Sides starting + * from the startSide + */ + //================================================================================ + + bool isClosedChainOfSides( _BlockSide* startSide, + map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides ) + { + // get start and end nodes + const SMDS_MeshNode *n1 = 0, *n2 = 0, *n; + for ( int y = 0; y < 2; ++y ) + for ( int x = 0; x < 2; ++x ) + { + n = startSide->getCornerNode(x,y); + if ( !corner2Sides.count( n )) continue; + if ( n1 ) + n2 = n; + else + n1 = n; + } + if ( !n2 ) return false; + + map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator + c2sides = corner2Sides.find( n1 ); + if ( c2sides == corner2Sides.end() ) return false; + + int nbChainLinks = 1; + n = n1; + _BlockSide* prevSide = startSide; + while ( n != n2 ) + { + // get the next side sharing n + list< _BlockSide* > & sides = c2sides->second; + _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() ); + if ( nextSide == prevSide ) return false; + + // find the next corner of the nextSide being in corner2Sides + n1 = n; + n = 0; + for ( int y = 0; y < 2 && !n; ++y ) + for ( int x = 0; x < 2; ++x ) + { + n = nextSide->getCornerNode(x,y); + c2sides = corner2Sides.find( n ); + if ( n == n1 || c2sides == corner2Sides.end() ) + n = 0; + else + break; + } + if ( !n ) return false; + + prevSide = nextSide; + nbChainLinks++; + } + + return ( n == n2 && nbChainLinks == NB_QUAD_SIDES ); + } + //================================================================================ /*! * \brief Try to find a block side adjacent to the given side by given edge */ //================================================================================ - _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide, - EQuadEdge sharedSideEdge1, - EQuadEdge sharedSideEdge2, - bool withGeometricAnalysis) + _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide, + EQuadEdge sharedSideEdge1, + EQuadEdge sharedSideEdge2, + bool withGeometricAnalysis, + set< _BlockSide* >& sidesAround) { _Block& block = _blocks.back(); _OrientedBlockSide& side1 = block._side[ startBlockSide ]; @@ -737,45 +792,75 @@ namespace 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 ) + if ( !withGeometricAnalysis ) { - _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); + sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() ); + return 0; } 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; + // Issue 0021529. There are at least 2 sides by each edge and + // position of block gravity center is undefined. + // Find a side starting from which we can walk around the startBlockSide + + // fill in corner2Sides + map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides; + for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt ) + { + _BlockSide* sideI = *sideIt; + corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI ); + corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI ); + corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI ); + corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI ); + } + // remove corners of startBlockSide from corner2Sides + set::iterator nIt = block._corners.begin(); + for ( ; nIt != block._corners.end(); ++nIt ) + corner2Sides.erase( *nIt ); + + // select a side + for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt ) + { + if ( isClosedChainOfSides( *sideIt, corner2Sides )) + { + foundSide = *sideIt; + break; + } + } + if ( !foundSide ) + return 0; } else { + // 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_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_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. * M_PI; // angle [0-2*PI] + angleOfSide.insert( make_pair( angle, sideI )); + _DUMP_(" "<< sideI << " - side dir (" + << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")" + << " angle " << angle); + } + 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 ] ) @@ -849,14 +934,12 @@ namespace row2.push_back( n1 = oppositeNode( quad, i1 )); } - _NodeDescriptor nDesc( row1[1] ); - if ( nDesc == _NodeDescriptor( row1[0] )) - return true; // row of 2 nodes + if ( isCornerNode( row1[1] )) + return true; // Find the rest nodes TIDSortedElemSet emptySet, avoidSet; - //while ( !isCornerNode( n2 )) - while ( nDesc == _NodeDescriptor( n2 )) + while ( !isCornerNode( n2 ) ) { avoidSet.clear(); avoidSet.insert( quad ); quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 ); @@ -1009,7 +1092,7 @@ bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* // ---------------------------- // Add internal nodes of a box // ---------------------------- - // projection points of internal nodes on box subshapes by which + // projection points of internal nodes on box sub-shapes by which // coordinates of internal nodes are computed vector pointOnShape( SMESH_Block::ID_Shell ); @@ -1218,4 +1301,3 @@ bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&) { return error("Algorithm can't work with geometrical shapes"); } -