From 9ed324256590c788619fc596b5c548799343e0db Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 24 Feb 2010 12:46:38 +0000 Subject: [PATCH] 0020611: [CEA] Algo lacking when editing MED meshing * Generate hexahedrons from 2D skin mesh: + StdMeshers_HexaFromSkin_3D.cxx --- src/StdMeshers/Makefile.am | 6 +- src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx | 1038 +++++++++++++++++ src/StdMeshers/StdMeshers_HexaFromSkin_3D.hxx | 57 + src/StdMeshers/StdMeshers_Hexa_3D.cxx | 22 +- src/StdMeshers/StdMeshers_Hexa_3D.hxx | 6 +- 5 files changed, 1123 insertions(+), 6 deletions(-) create mode 100644 src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx create mode 100644 src/StdMeshers/StdMeshers_HexaFromSkin_3D.hxx diff --git a/src/StdMeshers/Makefile.am b/src/StdMeshers/Makefile.am index c10d17e66..deeb6c2e7 100644 --- a/src/StdMeshers/Makefile.am +++ b/src/StdMeshers/Makefile.am @@ -73,7 +73,8 @@ salomeinclude_HEADERS = \ StdMeshers_CompositeHexa_3D.hxx \ StdMeshers_MaxLength.hxx \ StdMeshers_QuadrangleParams.hxx \ - StdMeshers_RadialQuadrangle_1D2D.hxx + StdMeshers_RadialQuadrangle_1D2D.hxx \ + StdMeshers_HexaFromSkin_3D.hxx # Libraries targets @@ -123,7 +124,8 @@ dist_libStdMeshers_la_SOURCES = \ StdMeshers_CompositeHexa_3D.cxx \ StdMeshers_MaxLength.cxx \ StdMeshers_QuadrangleParams.cxx \ - StdMeshers_RadialQuadrangle_1D2D.cxx + StdMeshers_RadialQuadrangle_1D2D.cxx \ + StdMeshers_HexaFromSkin_3D.cxx # additionnal information to compil and link file diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx new file mode 100644 index 000000000..464fc4c3e --- /dev/null +++ b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx @@ -0,0 +1,1038 @@ +// Copyright (C) 2007-2008 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 +// +// 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. +// +// 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 +// +// 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" +#include "SMDS_VolumeTool.hxx" +#include "SMESH_Block.hxx" +#include "SMESH_MesherHelper.hxx" + +#include "utilities.h" + +// Define error message +#ifdef _DEBUG_ +#define BAD_MESH_ERR \ + error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \ + __FILE__ ":" )<<__LINE__) +#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 + + + +namespace +{ + enum EBoxSides //!< sides of the block + { + B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED + }; + const char* SBoxSides[] = //!< names of block sides + { + "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED" + }; + enum EQuadEdge //!< edges of quadrangle side + { + Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED + }; + + //================================================================================ + /*! + * \brief return true if a node is at block corner + */ + //================================================================================ + + bool isCornerNode( const SMDS_MeshNode* n ) + { + return n && n->NbInverseElements( SMDSAbs_Face ) % 2; + } + + //================================================================================ + /*! + * \brief check element type + */ + //================================================================================ + + bool isQuadrangle(const SMDS_MeshElement* e) + { + return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 )); + } + + //================================================================================ + /*! + * \brief return opposite node of a quadrangle face + */ + //================================================================================ + + const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode) + { + return quad->GetNode( (iNode+2) % 4 ); + } + + //================================================================================ + /*! + * \brief Convertor of a pair of integers to a sole index + */ + struct _Indexer + { + int _xSize, _ySize; + _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {} + int size() const { return _xSize * _ySize; } + int operator()(int x, int y) const { return y * _xSize + x; } + }; + //================================================================================ + /*! + * \brief Oriented convertor of a pair of integers to a sole index + */ + class _OrientedIndexer : public _Indexer + { + public: + enum OriFlags //!< types of block side orientation + { + REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY + }; + _OrientedIndexer( const _Indexer& indexer, const int oriFlags ): + _Indexer( indexer._xSize, indexer._ySize ), + _xSize (indexer._xSize), _ySize(indexer._ySize), + _xRevFun((oriFlags & REV_X) ? & reverse : & lazy), + _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy), + _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy) + { + (*_swapFun)( _xSize, _ySize ); + } + //!< Return index by XY + int operator()(int x, int y) const + { + (*_xRevFun)( x, const_cast( _xSize )); + (*_yRevFun)( y, const_cast( _ySize )); + (*_swapFun)( x, y ); + return _Indexer::operator()(x,y); + } + //!< Return index for a corner + int corner(bool xMax, bool yMax) const + { + int x = xMax, y = yMax, size = 2; + (*_xRevFun)( x, size ); + (*_yRevFun)( y, size ); + (*_swapFun)( x, y ); + return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0); + } + int xSize() const { return _xSize; } + int ySize() const { return _ySize; } + private: + _Indexer _indexer; + int _xSize, _ySize; + + typedef void (*TFun)(int& x, int& y); + TFun _xRevFun, _yRevFun, _swapFun; + + 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); } + }; + //================================================================================ + /*! + * \brief Structure corresponding to the meshed side of block + */ + struct _BlockSide + { + vector _grid; + _Indexer _index; + int _nbBlocksExpected; + int _nbBlocksFound; + +#ifdef _DEBUG_ +#define _grid_access_(args) _grid.at( args ) +#else +#define _grid_access_(args) _grid[ args ] +#endif + //!< Return node at XY + const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); } + //!< Set node at XY + void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; } + //!< Return a corner node + const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const + { + 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 + 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 )); } + //!< Return gravity center of the four corners and the middle node + gp_XYZ getGC() const + { + gp_XYZ xyz = + getXYZ( 0, 0 ) + + getXYZ( _index._xSize-1, 0 ) + + getXYZ( 0, _index._ySize-1 ) + + getXYZ( _index._xSize-1, _index._ySize-1 ) + + getXYZ( _index._xSize/2, _index._ySize/2 ); + return xyz / 5; + } + //!< Return number of mesh faces + int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); } + }; + //================================================================================ + /*! + * \brief _BlockSide with changed orientation + */ + struct _OrientedBlockSide + { + _BlockSide* _side; + _OrientedIndexer _index; + + _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ): + _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {} + //!< return coordinates by XY + gp_XYZ xyz(int x, int y) const + { + return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) ); + } + //!< safely return a node by XY + const SMDS_MeshNode* node(int x, int y) const + { + int i = _index( x, y ); + return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i]; + } + //!< Return a corner node + const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const + { + return _side->_grid_access_( _index.corner( isXMax, isYMax )); + } + //!< return its size in nodes + int getHoriSize() const { return _index.xSize(); } + int getVertSize() const { return _index.ySize(); } + //!< True if _side has been initialized + operator bool() const { return _side; } + //! Direct access to _side + const _BlockSide* operator->() const { return _side; } + _BlockSide* operator->() { return _side; } + }; + //================================================================================ + /*! + * \brief Meshed skin of block + */ + struct _Block + { + _OrientedBlockSide _side[6]; // 6 sides of a sub-block + + const _OrientedBlockSide& getSide(int i) const { return _side[i]; } + }; + //================================================================================ + /*! + * \brief Skin mesh possibly containing several meshed blocks + */ + class _Skin + { + public: + + int findBlocks(SMESH_Mesh& mesh); + //!< return i-th block + const _Block& getBlock(int i) const { return _blocks[i]; } + //!< return error description + const SMESH_Comment& error() const { return _error; } + + private: + bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad); + bool fillRowsUntilCorner(const SMDS_MeshElement* quad, + const SMDS_MeshNode* n1, + const SMDS_MeshNode* n2, + vector& verRow1, + vector& verRow2, + bool alongN1N2 ); + _OrientedBlockSide findBlockSide( EBoxSides startBlockSide, + EQuadEdge sharedSideEdge1, + EQuadEdge sharedSideEdge2); + //!< update own data and data of the side bound to block + void setSideBoundToBlock( _BlockSide& side ) + { + side._nbBlocksFound++; + if ( side.isBound() ) + { + _corner2sides[ side.getCornerNode(0,0) ].erase( &side ); + _corner2sides[ side.getCornerNode(1,0) ].erase( &side ); + _corner2sides[ side.getCornerNode(0,1) ].erase( &side ); + _corner2sides[ side.getCornerNode(1,1) ].erase( &side ); + } + } + //!< store reason of error + int error(const SMESH_Comment& reason) { _error = reason; return 0; } + + SMESH_Comment _error; + + list< _BlockSide > _allSides; + vector< _Block > _blocks; + + map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides; + }; + + //================================================================================ + /*! + * \brief Find and return number of submeshes corresponding to blocks + */ + //================================================================================ + + int _Skin::findBlocks(SMESH_Mesh& mesh) + { + SMESHDS_Mesh* meshDS = mesh.GetMeshDS(); + + // Find a node at any block corner + + SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true); + if ( !nIt->more() ) return error("Empty mesh"); + + const SMDS_MeshNode* nCorner = 0; + while ( nIt->more() ) + { + nCorner = nIt->next(); + if ( isCornerNode( nCorner )) + break; + else + nCorner = 0; + } + if ( !nCorner ) + return BAD_MESH_ERR; + + // -------------------------------------------------------------------- + // Find all block sides starting from mesh faces sharing the corner node + // -------------------------------------------------------------------- + + int nbFacesOnSides = 0; + TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's + list< const SMDS_MeshNode* > corners( 1, nCorner ); + list< const SMDS_MeshNode* >::iterator corner = corners.begin(); + while ( corner != corners.end() ) + { + SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face ); + while ( faceIt->more() ) + { + const SMDS_MeshElement* face = faceIt->next(); + if ( !cornerFaces.insert( face ).second ) + continue; // already loaded block side + + if ( !isQuadrangle( face )) + return error("Non-quadrangle elements in the input mesh"); + + if ( _allSides.empty() || !_allSides.back()._grid.empty() ) + _allSides.push_back( _BlockSide() ); + + _BlockSide& side = _allSides.back(); + if ( !fillSide( side, face ) ) + { + if ( !_error.empty() ) + return false; + } + else + { + for ( int isXMax = 0; isXMax < 2; ++isXMax ) + 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 )); + } + nbFacesOnSides += side.getNbFaces(); + } + } + ++corner; + + // find block sides of other domains if any + if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() ) + { + while ( nIt->more() ) + { + nCorner = nIt->next(); + if ( isCornerNode( nCorner )) + corner = corners.insert( corner, nCorner ); + } + nbFacesOnSides = mesh.NbQuadrangles(); + } + } + + if ( _allSides.empty() ) + return BAD_MESH_ERR; + if ( _allSides.back()._grid.empty() ) + _allSides.pop_back(); + + // --------------------------- + // 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 ) + { + _BlockSide& side = *sideIt; + bool isSharedSide = true; + for ( int isXMax = 0; isXMax < 2; ++isXMax ) + for ( int isYMax = 0; isYMax < 2; ++isYMax ) + if ( _corner2sides[ side.getCornerNode(isXMax,isYMax) ].size() < 5 ) + isSharedSide = false; + side._nbBlocksFound = 0; + side._nbBlocksExpected = isSharedSide ? 2 : 1; + nbBlockSides += side._nbBlocksExpected; + } + + // find sides of each block + int nbBlocks = 0; + while ( nbBlockSides >= 6 ) + { + // get any side not bound to all blocks to belongs to + sideIt = _allSides.begin(); + while ( sideIt != _allSides.end() && sideIt->isBound()) + ++sideIt; + + // start searching for block sides from the got side + bool ok = true; + if ( _blocks.empty() || _blocks.back()._side[B_FRONT] ) + _blocks.resize( _blocks.size() + 1 ); + + _Block& block = _blocks.back(); + block._side[B_FRONT] = &(*sideIt); + setSideBoundToBlock( *sideIt ); + 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 )); + + // count the found sides + _DUMP_(endl); + for (int i = 0; i < B_UNDEFINED; ++i ) + { + _DUMP_("** Block side "<< SBoxSides[i]); + if ( block._side[ i ] ) + { + if ( ok && i != B_FRONT) + { + 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)); + } + else + { + _DUMP_("Not found"<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); + + // find out size of block side + vector horRow1, horRow2, verRow1, verRow2; + if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) || + !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false )) + return false; + nbX = horRow1.size(), nbY = verRow1.size(); + + // store found nodes + side._index._xSize = horRow1.size(); + side._index._ySize = verRow1.size(); + side._grid.resize( side._index.size(), NULL ); + + for ( x = 0; x < horRow1.size(); ++x ) + { + side.setNode( x, 0, horRow1[x] ); + side.setNode( x, 1, horRow2[x] ); + } + for ( y = 0; y < verRow1.size(); ++y ) + { + side.setNode( 0, y, verRow1[y] ); + side.setNode( 1, y, verRow2[y] ); + } + } + // Find the rest nodes + + y = 1; // y of the row to fill + TIDSortedElemSet emptySet, avoidSet; + while ( ++y < nbY ) + { + // get next firstQuad in the next row of quadrangles + // + // n2up + // o---o <- y row + // | | + // o---o o o o o <- found nodes + //n1down n2down + // + int i1down, i2down, i2up; + const SMDS_MeshNode* n1down = side.getNode( 0, y-1 ); + const SMDS_MeshNode* n2down = side.getNode( 1, y-1 ); + avoidSet.clear(); avoidSet.insert( firstQuad ); + firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet, + &i1down, &i2down); + if ( !isQuadrangle( firstQuad )) + return BAD_MESH_ERR; + + const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down ); + avoidSet.clear(); avoidSet.insert( firstQuad ); + + // find the rest nodes in the y-th row by faces in the row + + x = 1; + while ( ++x < nbX ) + { + const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet, + avoidSet, &i2up, &i2down); + if ( !isQuadrangle( quad )) + return BAD_MESH_ERR; + + n2up = oppositeNode( quad, i2down ); + n2down = oppositeNode( quad, i2up ); + avoidSet.clear(); avoidSet.insert( quad ); + + side.setNode( x, y, n2up ); + } + } + + return true; + } + + //================================================================================ + /*! + * \brief Try to find a block side adjacent to the given side by given edge + */ + //================================================================================ + + _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide, + EQuadEdge sharedSideEdge1, + EQuadEdge sharedSideEdge2) + { + _Block& block = _blocks.back(); + _OrientedBlockSide& side1 = block._side[ startBlockSide ]; + + // get corner nodes of the given block edge + bool xMax1=0, yMax1=0, xMax2=1, yMax2=1; + switch( sharedSideEdge1 ) + { + case Q_BOTTOM: yMax2 = 0; break; + case Q_RIGHT: xMax1 = 1; break; + case Q_TOP: yMax1 = 1; break; + case Q_LEFT: xMax2 = 0; break; + default: + error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__); + return 0; + } + const SMDS_MeshNode* n1 = side1.cornerNode( xMax1, yMax1); + const SMDS_MeshNode* n2 = side1.cornerNode( xMax2, yMax2); + + set< _BlockSide* >& sidesWithN1 = _corner2sides[ n1 ]; + set< _BlockSide* >& sidesWithN2 = _corner2sides[ n2 ]; + if ( sidesWithN1.empty() || sidesWithN2.empty() ) + { + _DUMP_("no sides by nodes "<< n1->GetID() << " " << n2->GetID() ); + return 0; + } + + // find all sides sharing both nodes n1 and n2 + set< _BlockSide* > sidesOnEdge; + set< _BlockSide* >::iterator sideIt; + std::set_intersection( sidesWithN1.begin(), sidesWithN1.end(), + sidesWithN2.begin(), sidesWithN2.end(), + inserter( sidesOnEdge, sidesOnEdge.begin())); + + // exclude loaded sides of block from sidesOnEdge + int nbLoadedSides = 0; + for (int i = 0; i < B_UNDEFINED; ++i ) + { + if ( block._side[ i ] ) + { + nbLoadedSides++; + sidesOnEdge.erase( block._side[ i ]._side ); + } + } + int nbSidesOnEdge = sidesOnEdge.size(); + _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge); + if ( nbSidesOnEdge == 0 ) + return 0; + + _BlockSide* foundSide = 0; + if ( nbSidesOnEdge == 1 || nbSidesOnEdge == 2 && nbLoadedSides == 1 ) + { + foundSide = *sidesOnEdge.begin(); + } + else + { + // Select one of found sides most close to startBlockSide + + // gravity center of already loaded block sides + gp_XYZ gc; + for (int i = 0; i < B_UNDEFINED; ++i ) + if ( block._side[ i ] ) + { + gc += block._side[ i ]._side->getGC(); + nbLoadedSides++; + } + gc /= nbLoadedSides; + + gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z()); + gp_Vec p2p1( 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 ); + + 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 ); + double angle = side1Dir.AngleWithRef( sideIDir, p2p1 ); + if ( angle < 0 ) angle += 2 * PI; + angleOfSide.insert( make_pair( angle, sideI )); + } + gp_Vec gcDir( p1, gc ); + double gcAngle = side1Dir.AngleWithRef( gcDir, p2p1 ); + foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second; + } + + // Orient the found side correctly + + // corners of found side corresponding to nodes n1 and n2 + xMax1= yMax1 = 0; xMax2 = yMax2 = 1; + switch( sharedSideEdge2 ) + { + case Q_BOTTOM: yMax2 = 0; break; + case Q_RIGHT: xMax1 = 1; break; + case Q_TOP: yMax1 = 1; break; + case Q_LEFT: xMax2 = 0; break; + default: + error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__); + return 0; + } + for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori ) + { + _OrientedBlockSide orientedSide( foundSide, ori ); + const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1); + const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2); + if ( n1 == n12 && n2 == n22 ) + return orientedSide; + } + error(SMESH_Comment("Failed to orient a block side found by edge ")<& row1, + vector& row2, + const bool alongN1N2 ) + { + const SMDS_MeshNode* corner1 = n1; + + // Store nodes of quad in the rows and find new n1 and n2 to get + // the next face so that new n2 is on block edge + int i1 = quad->GetNodeIndex( n1 ); + int i2 = quad->GetNodeIndex( n2 ); + row1.clear(); row2.clear(); + row1.push_back( n1 ); + if ( alongN1N2 ) + { + row1.push_back( n2 ); + row2.push_back( oppositeNode( quad, i2 )); + row2.push_back( n1 = oppositeNode( quad, i1 )); + } + else + { + row2.push_back( n2 ); + row1.push_back( n2 = oppositeNode( quad, i2 )); + row2.push_back( n1 = oppositeNode( quad, i1 )); + } + + // Find the rest nodes + TIDSortedElemSet emptySet, avoidSet; + while ( !isCornerNode( n2 )) + { + avoidSet.clear(); avoidSet.insert( quad ); + quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 ); + if ( !isQuadrangle( quad )) + return BAD_MESH_ERR; + + row1.push_back( n2 = oppositeNode( quad, i1 )); + row2.push_back( n1 = oppositeNode( quad, i2 )); + } + return n1 != corner1; + } + + //================================================================================ + /*! + * \brief Return a corner face by a corner node + */ + //================================================================================ + + const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const + { + int x, y, isXMax, isYMax, found = 0; + for ( isXMax = 0; isXMax < 2; ++isXMax ) + { + for ( isYMax = 0; isYMax < 2; ++isYMax ) + { + x = isXMax ? _index._xSize-1 : 0; + y = isYMax ? _index._ySize-1 : 0; + found = ( getNode(x,y) == cornerNode ); + if ( found ) break; + } + if ( found ) break; + } + if ( !found ) return 0; + int dx = isXMax ? -1 : +1; + int dy = isYMax ? -1 : +1; + const SMDS_MeshNode* n1 = getNode(x,y); + const SMDS_MeshNode* n2 = getNode(x+dx,y); + const SMDS_MeshNode* n3 = getNode(x,y+dy); + const SMDS_MeshNode* n4 = getNode(x+dx,y+dy); + return SMDS_Mesh::FindFace(n1, n2, n3, n4 ); + } + +} + +//======================================================================= +//function : StdMeshers_HexaFromSkin_3D +//purpose : +//======================================================================= + +StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen) + :SMESH_3D_Algo(hypId, studyId, gen) +{ + MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D"); + _name = "HexaFromSkin_3D"; +} + +StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D() +{ + MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D"); +} + +//================================================================================ +/*! + * \brief Main method, which generates hexaheda + */ +//================================================================================ + +bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper) +{ + _Skin skin; + int nbBlocks = skin.findBlocks(aMesh); + if ( nbBlocks == 0 ) + return error( skin.error()); + + vector< vector< const SMDS_MeshNode* > > columns; + int x, xSize, y, ySize, z, zSize; + _Indexer colIndex; + + for ( int i = 0; i < nbBlocks; ++i ) + { + const _Block& block = skin.getBlock( i ); + + // ------------------------------------------ + // Fill columns of nodes with existing nodes + // ------------------------------------------ + + xSize = block.getSide(B_BOTTOM).getHoriSize(); + ySize = block.getSide(B_BOTTOM).getVertSize(); + zSize = block.getSide(B_FRONT ).getVertSize(); + int X = xSize - 1, Y = ySize - 1, Z = zSize - 1; + colIndex = _Indexer( xSize, ySize ); + columns.resize( colIndex.size() ); + + // fill node columns by front and back box sides + for ( x = 0; x < xSize; ++x ) { + vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )]; + vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )]; + column0.resize( zSize ); + column1.resize( zSize ); + for ( z = 0; z < zSize; ++z ) { + column0[ z ] = block.getSide(B_FRONT).node( x, z ); + column1[ z ] = block.getSide(B_BACK) .node( x, z ); + } + } + // fill node columns by left and right box sides + for ( y = 1; y < ySize-1; ++y ) { + vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )]; + vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )]; + column0.resize( zSize ); + column1.resize( zSize ); + for ( z = 0; z < zSize; ++z ) { + column0[ z ] = block.getSide(B_LEFT) .node( y, z ); + column1[ z ] = block.getSide(B_RIGHT).node( y, z ); + } + } + // get nodes from top and bottom box sides + for ( x = 1; x < xSize-1; ++x ) { + for ( y = 1; y < ySize-1; ++y ) { + vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )]; + column.resize( zSize ); + column.front() = block.getSide(B_BOTTOM).node( x, y ); + column.back() = block.getSide(B_TOP) .node( x, y ); + } + } + + // ---------------------------- + // Add internal nodes of a box + // ---------------------------- + // projection points of internal nodes on box subshapes by which + // coordinates of internal nodes are computed + vector pointOnShape( SMESH_Block::ID_Shell ); + + // projections on vertices are constant + pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 ); + pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 ); + pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y ); + pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y ); + pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 ); + pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 ); + pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y ); + pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y ); + + for ( x = 1; x < xSize-1; ++x ) + { + gp_XYZ params; // normalized parameters of internal node within a unit box + params.SetCoord( 1, x / double(X) ); + for ( y = 1; y < ySize-1; ++y ) + { + params.SetCoord( 2, y / double(Y) ); + // column to fill during z loop + vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )]; + // projections on horizontal edges + pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 ); + pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y ); + pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y ); + pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y ); + pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 ); + pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y ); + pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y ); + pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y ); + // projections on horizontal sides + pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y ); + pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y ); + for ( z = 1; z < zSize-1; ++z ) // z loop + { + params.SetCoord( 3, z / double(Z) ); + // projections on vertical edges + pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z ); + pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z ); + pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z ); + pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z ); + // projections on vertical sides + pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z ); + pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z ); + pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z ); + pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z ); + + // compute internal node coordinates + gp_XYZ coords; + SMESH_Block::ShellPoint( params, pointOnShape, coords ); + column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() ); + +#ifdef DEB_GRID + // debug + //cout << "----------------------------------------------------------------------"<& col00 = columns[ colIndex( x, y )]; + vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )]; + vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )]; + vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )]; + // bottom face normal of a hexa mush point outside the volume + if ( isForw ) + for ( z = 0; z < zSize-1; ++z ) + aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z], + col00[z+1], col01[z+1], col11[z+1], col10[z+1]); + else + for ( z = 0; z < zSize-1; ++z ) + aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z], + col00[z+1], col10[z+1], col11[z+1], col01[z+1]); + } + } + } // loop on blocks + + return true; +} + +//================================================================================ +/*! + * \brief Evaluate nb of hexa + */ +//================================================================================ + +bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap) +{ + _Skin skin; + int nbBlocks = skin.findBlocks(aMesh); + if ( nbBlocks == 0 ) + return error( skin.error()); + + bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC ); + + int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa; + vector& nbByType = aResMap[ aMesh.GetSubMesh( aShape )]; + if ( entity >= nbByType.size() ) + nbByType.resize( SMDSEntity_Last, 0 ); + + for ( int i = 0; i < nbBlocks; ++i ) + { + const _Block& block = skin.getBlock( i ); + + int nbX = block.getSide(B_BOTTOM).getHoriSize(); + int nbY = block.getSide(B_BOTTOM).getVertSize(); + int nbZ = block.getSide(B_FRONT ).getVertSize(); + + int nbHexa = (nbX-1) * (nbY-1) * (nbZ-1); + int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2); + if ( secondOrder ) + nbNodes += + (nbX-2) * (nbY-2) * (nbZ-1) + + (nbX-2) * (nbY-1) * (nbZ-2) + + (nbX-1) * (nbY-2) * (nbZ-2); + + + nbByType[ entity ] += nbHexa; + nbByType[ SMDSEntity_Node ] += nbNodes; + } + + return true; +} + +//================================================================================ +/*! + * \brief Abstract method must be defined but does nothing + */ +//================================================================================ + +bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&, + Hypothesis_Status& aStatus) +{ + aStatus = SMESH_Hypothesis::HYP_OK; + return true; +} + +//================================================================================ +/*! + * \brief Abstract method must be defined but just reports an error as this + * algo is not intended to work with shapes + */ +//================================================================================ + +bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&) +{ + return error("Algorithm can't work with geometrical shapes"); +} + diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.hxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.hxx new file mode 100644 index 000000000..32126da6b --- /dev/null +++ b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.hxx @@ -0,0 +1,57 @@ +// Copyright (C) 2007-2008 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 +// +// 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. +// +// 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 +// +// File : StdMeshers_HexaFromSkin_3D.hxx +// Created : Wed Jan 27 12:23:21 2010 +// Author : Edward AGAPOV (eap) + + +#ifndef __StdMeshers_HexaFromSkin_3D_HXX__ +#define __StdMeshers_HexaFromSkin_3D_HXX__ + +#include "SMESH_StdMeshers.hxx" +#include "SMESH_3D_Algo.hxx" + +/*! + * \brief Alorithm generating hexahedral mesh from 2D skin of block + */ + +class STDMESHERS_EXPORT StdMeshers_HexaFromSkin_3D : public SMESH_3D_Algo +{ +public: + StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen); + virtual ~StdMeshers_HexaFromSkin_3D(); + + virtual bool Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper); + + virtual bool CheckHypothesis(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus); + + virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + + virtual bool Evaluate(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap); + +}; + +#endif diff --git a/src/StdMeshers/StdMeshers_Hexa_3D.cxx b/src/StdMeshers/StdMeshers_Hexa_3D.cxx index 42b5488ec..0e68f0214 100644 --- a/src/StdMeshers/StdMeshers_Hexa_3D.cxx +++ b/src/StdMeshers/StdMeshers_Hexa_3D.cxx @@ -26,8 +26,10 @@ // Module : SMESH // #include "StdMeshers_Hexa_3D.hxx" + #include "StdMeshers_CompositeHexa_3D.hxx" #include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_HexaFromSkin_3D.hxx" #include "StdMeshers_Penta_3D.hxx" #include "StdMeshers_Prism_3D.hxx" #include "StdMeshers_Quadrangle_2D.hxx" @@ -78,6 +80,7 @@ StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen) MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D"); _name = "Hexa_3D"; _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID); // 1 bit /shape type + _requireShape = false; } //============================================================================= @@ -862,10 +865,27 @@ bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh, return true; } +//================================================================================ +/*! + * \brief Computes hexahedral mesh from 2D mesh of block + */ +//================================================================================ + +bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper) +{ + static StdMeshers_HexaFromSkin_3D * algo = 0; + if ( !algo ) { + SMESH_Gen* gen = aMesh.GetGen(); + algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen ); + } + algo->InitComputeError(); + algo->Compute( aMesh, aHelper ); + return error( algo->GetComputeError()); +} //============================================================================= /*! - * + * */ //============================================================================= diff --git a/src/StdMeshers/StdMeshers_Hexa_3D.hxx b/src/StdMeshers/StdMeshers_Hexa_3D.hxx index c892ac729..1e4c9a2fe 100644 --- a/src/StdMeshers/StdMeshers_Hexa_3D.hxx +++ b/src/StdMeshers/StdMeshers_Hexa_3D.hxx @@ -73,9 +73,9 @@ public: const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus); - virtual bool Compute(SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape) - /*throw (SALOME_Exception)*/; + virtual bool Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape); + + virtual bool Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper); virtual bool Evaluate(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape, MapShapeNbElems& aResMap);