src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : StdMeshers_HexaFromSkin_3D.cxx
  21 // Created   : Wed Jan 27 12:28:07 2010
  22 // Author    : Edward AGAPOV (eap)
  23 //
  24 #include "StdMeshers_HexaFromSkin_3D.hxx"
  25
  26 #include "SMDS_VolumeOfNodes.hxx"
  27 #include "SMDS_VolumeTool.hxx"
  28 #include "SMESH_Block.hxx"
  29 #include "SMESH_MesherHelper.hxx"
  30
  31 #include "utilities.h"
  32
  33 // Define error message
  34 #ifdef _DEBUG_
  35 #define BAD_MESH_ERR \
  36   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
  37                       __FILE__ ":" )<<__LINE__)
  38 #else
  39 #define BAD_MESH_ERR \
  40   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
  41 #endif
  42
  43 // Debug output
  44 #define _DUMP_(msg) // cout << msg << endl
  45
  46
  47
  48 namespace
  49 {
  50   enum EBoxSides //!< sides of the block
  51     {
  52       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED
  53     };
  54   const char* SBoxSides[] = //!< names of block sides
  55     {
  56       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
  57     };
  58   enum EQuadEdge //!< edges of quadrangle side
  59     {
  60       Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED
  61     };
  62
  63   //================================================================================
  64   /*!
  65    * \brief return true if a node is at block corner
  66    */
  67   //================================================================================
  68
  69   bool isCornerNode( const SMDS_MeshNode* n )
  70   {
  71     return n && n->NbInverseElements( SMDSAbs_Face ) % 2;
  72   }
  73
  74   //================================================================================
  75   /*!
  76    * \brief check element type
  77    */
  78   //================================================================================
  79
  80   bool isQuadrangle(const SMDS_MeshElement* e)
  81   {
  82     return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 ));
  83   }
  84
  85   //================================================================================
  86   /*!
  87    * \brief return opposite node of a quadrangle face
  88    */
  89   //================================================================================
  90
  91   const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
  92   {
  93     return quad->GetNode( (iNode+2) % 4 );
  94   }
  95
  96   //================================================================================
  97   /*!
  98    * \brief Convertor of a pair of integers to a sole index
  99    */
 100   struct _Indexer
 101   {
 102     int _xSize, _ySize;
 103     _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
 104     int size() const { return _xSize * _ySize; }
 105     int operator()(int x, int y) const { return y * _xSize + x; }
 106   };
 107   //================================================================================
 108   /*!
 109    * \brief Oriented convertor of a pair of integers to a sole index
 110    */
 111   class _OrientedIndexer : public _Indexer
 112   {
 113   public:
 114     enum OriFlags //!< types of block side orientation
 115       {
 116         REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
 117       };
 118     _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
 119       _Indexer( indexer._xSize, indexer._ySize ),
 120       _xSize (indexer._xSize), _ySize(indexer._ySize),
 121       _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
 122       _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
 123       _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
 124     {
 125       (*_swapFun)( _xSize, _ySize );
 126     }
 127     //!< Return index by XY
 128     int operator()(int x, int y) const
 129     {
 130       (*_xRevFun)( x, const_cast<int&>( _xSize ));
 131       (*_yRevFun)( y, const_cast<int&>( _ySize ));
 132       (*_swapFun)( x, y );
 133       return _Indexer::operator()(x,y);
 134     }
 135     //!< Return index for a corner
 136     int corner(bool xMax, bool yMax) const
 137     {
 138       int x = xMax, y = yMax, size = 2;
 139       (*_xRevFun)( x, size );
 140       (*_yRevFun)( y, size );
 141       (*_swapFun)( x, y );
 142       return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
 143     }
 144     int xSize() const { return _xSize; }
 145     int ySize() const { return _ySize; }
 146   private:
 147     _Indexer _indexer;
 148     int _xSize, _ySize;
 149
 150     typedef void (*TFun)(int& x, int& y);
 151     TFun _xRevFun, _yRevFun, _swapFun;
 152
 153     static void lazy(int&, int&) {}
 154     static void reverse(int& x, int& size) { x = size - x - 1; }
 155     static void swap(int& x, int& y) { std::swap(x,y); }
 156   };
 157   //================================================================================
 158   /*!
 159    * \brief Structure corresponding to the meshed side of block
 160    */
 161   struct _BlockSide
 162   {
 163     vector<const SMDS_MeshNode*> _grid;
 164     _Indexer                     _index;
 165     int                          _nbBlocksExpected;
 166     int                          _nbBlocksFound;
 167
 168 #ifdef _DEBUG_
 169 #define _grid_access_(args) _grid.at( args )
 170 #else
 171 #define _grid_access_(args) _grid[ args ]
 172 #endif
 173     //!< Return node at XY
 174     const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); }
 175     //!< Set node at XY
 176     void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; }
 177     //!< Return a corner node
 178     const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
 179     {
 180       return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
 181     }
 182     const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
 183     //!< True if all blocks this side belongs to have beem found
 184     bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
 185     //!< Return coordinates of node at XY
 186     gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); }
 187     //!< Return gravity center of the four corners and the middle node
 188     gp_XYZ getGC() const
 189     {
 190       gp_XYZ xyz =
 191         getXYZ( 0, 0 ) +
 192         getXYZ( _index._xSize-1, 0 ) +
 193         getXYZ( 0, _index._ySize-1 ) +
 194         getXYZ( _index._xSize-1, _index._ySize-1 ) +
 195         getXYZ( _index._xSize/2, _index._ySize/2 );
 196       return xyz / 5;
 197     }
 198     //!< Return number of mesh faces
 199     int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
 200   };
 201   //================================================================================
 202   /*!
 203    * \brief _BlockSide with changed orientation
 204    */
 205   struct _OrientedBlockSide
 206   {
 207     _BlockSide*       _side;
 208     _OrientedIndexer  _index;
 209
 210     _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
 211       _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
 212     //!< return coordinates by XY
 213     gp_XYZ xyz(int x, int y) const
 214     {
 215       return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) );
 216     }
 217     //!< safely return a node by XY
 218     const SMDS_MeshNode* node(int x, int y) const
 219     {
 220       int i = _index( x, y );
 221       return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
 222     }
 223     //!< Return a corner node
 224     const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
 225     {
 226       return _side->_grid_access_( _index.corner( isXMax, isYMax ));
 227     }
 228     //!< return its size in nodes
 229     int getHoriSize() const { return _index.xSize(); }
 230     int getVertSize() const  { return _index.ySize(); }
 231     //!< True if _side has been initialized
 232     operator bool() const { return _side; }
 233     //! Direct access to _side
 234     const _BlockSide* operator->() const { return _side; }
 235     _BlockSide* operator->() { return _side; }
 236   };
 237   //================================================================================
 238   /*!
 239    * \brief Meshed skin of block
 240    */
 241   struct _Block
 242   {
 243     _OrientedBlockSide _side[6]; // 6 sides of a sub-block
 244
 245     const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
 246   };
 247   //================================================================================
 248   /*!
 249    * \brief Skin mesh possibly containing several meshed blocks
 250    */
 251   class _Skin
 252   {
 253   public:
 254
 255     int findBlocks(SMESH_Mesh& mesh);
 256     //!< return i-th block
 257     const _Block& getBlock(int i) const { return _blocks[i]; }
 258     //!< return error description
 259     const SMESH_Comment& error() const { return _error; }
 260
 261   private:
 262     bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad);
 263     bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
 264                              const SMDS_MeshNode*    n1,
 265                              const SMDS_MeshNode*    n2,
 266                              vector<const SMDS_MeshNode*>& verRow1,
 267                              vector<const SMDS_MeshNode*>& verRow2,
 268                              bool alongN1N2 );
 269     _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
 270                                       EQuadEdge sharedSideEdge1,
 271                                       EQuadEdge sharedSideEdge2);
 272     //!< update own data and data of the side bound to block
 273     void setSideBoundToBlock( _BlockSide& side )
 274     {
 275       side._nbBlocksFound++;
 276       if ( side.isBound() )
 277       {
 278         _corner2sides[ side.getCornerNode(0,0) ].erase( &side );
 279         _corner2sides[ side.getCornerNode(1,0) ].erase( &side );
 280         _corner2sides[ side.getCornerNode(0,1) ].erase( &side );
 281         _corner2sides[ side.getCornerNode(1,1) ].erase( &side );
 282       }
 283     }
 284     //!< store reason of error
 285     int error(const SMESH_Comment& reason) { _error = reason; return 0; }
 286
 287     SMESH_Comment      _error;
 288
 289     list< _BlockSide > _allSides;
 290     vector< _Block >   _blocks;
 291
 292     map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
 293   };
 294
 295   //================================================================================
 296   /*!
 297    * \brief Find and return number of submeshes corresponding to blocks
 298    */
 299   //================================================================================
 300
 301   int _Skin::findBlocks(SMESH_Mesh& mesh)
 302   {
 303     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 304
 305     // Find a node at any block corner
 306
 307     SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
 308     if ( !nIt->more() ) return error("Empty mesh");
 309
 310     const SMDS_MeshNode* nCorner = 0;
 311     while ( nIt->more() )
 312     {
 313       nCorner = nIt->next();
 314       if ( isCornerNode( nCorner ))
 315         break;
 316       else
 317         nCorner = 0;
 318     }
 319     if ( !nCorner )
 320       return BAD_MESH_ERR;
 321
 322     // --------------------------------------------------------------------
 323     // Find all block sides starting from mesh faces sharing the corner node
 324     // --------------------------------------------------------------------
 325
 326     int nbFacesOnSides = 0;
 327     TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
 328     list< const SMDS_MeshNode* > corners( 1, nCorner );
 329     list< const SMDS_MeshNode* >::iterator corner = corners.begin();
 330     while ( corner != corners.end() )
 331     {
 332       SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
 333       while ( faceIt->more() )
 334       {
 335         const SMDS_MeshElement* face = faceIt->next();
 336         if ( !cornerFaces.insert( face ).second )
 337           continue; // already loaded block side
 338
 339         if ( !isQuadrangle( face ))
 340           return error("Non-quadrangle elements in the input mesh");
 341
 342         if ( _allSides.empty() || !_allSides.back()._grid.empty() )
 343           _allSides.push_back( _BlockSide() );
 344
 345         _BlockSide& side = _allSides.back();
 346         if ( !fillSide( side, face ) )
 347         {
 348           if ( !_error.empty() )
 349             return false;
 350         }
 351         else
 352         {
 353           for ( int isXMax = 0; isXMax < 2; ++isXMax )
 354             for ( int isYMax = 0; isYMax < 2; ++isYMax )
 355             {
 356               const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
 357               if ( !isCornerNode( nCorner ))
 358                 return BAD_MESH_ERR;
 359               _corner2sides[ nCorner ].insert( &side );
 360               corners.push_back( nCorner );
 361               cornerFaces.insert( side.getCornerFace( nCorner ));
 362             }
 363           nbFacesOnSides += side.getNbFaces();
 364         }
 365       }
 366       ++corner;
 367
 368       // find block sides of other domains if any
 369       if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
 370       {
 371         while ( nIt->more() )
 372         {
 373           nCorner = nIt->next();
 374           if ( isCornerNode( nCorner ))
 375             corner = corners.insert( corner, nCorner );
 376         }
 377         nbFacesOnSides = mesh.NbQuadrangles();
 378       }
 379     }
 380
 381     if ( _allSides.empty() )
 382       return BAD_MESH_ERR;
 383     if ( _allSides.back()._grid.empty() )
 384       _allSides.pop_back();
 385
 386     // ---------------------------
 387     // Organize sides into blocks
 388     // ---------------------------
 389
 390     // analyse sharing of sides by blocks
 391     int nbBlockSides = 0; // nb of block sides taking into account their sharing
 392     list < _BlockSide >::iterator sideIt = _allSides.begin();
 393     for ( ; sideIt != _allSides.end(); ++sideIt )
 394     {
 395       _BlockSide& side = *sideIt;
 396       bool isSharedSide = true;
 397       for ( int isXMax = 0; isXMax < 2; ++isXMax )
 398         for ( int isYMax = 0; isYMax < 2; ++isYMax )
 399           if ( _corner2sides[ side.getCornerNode(isXMax,isYMax) ].size() < 5 )
 400             isSharedSide = false;
 401       side._nbBlocksFound = 0;
 402       side._nbBlocksExpected = isSharedSide ? 2 : 1;
 403       nbBlockSides += side._nbBlocksExpected;
 404     }
 405
 406     // find sides of each block
 407     int nbBlocks = 0;
 408     while ( nbBlockSides >= 6 )
 409     {
 410       // get any side not bound to all blocks it belongs to
 411       sideIt = _allSides.begin();
 412       while ( sideIt != _allSides.end() && sideIt->isBound())
 413         ++sideIt;
 414
 415       // start searching for block sides from the got side
 416       bool ok = true;
 417       if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
 418         _blocks.resize( _blocks.size() + 1 );
 419
 420       _Block& block = _blocks.back();
 421       block._side[B_FRONT] = &(*sideIt);
 422       setSideBoundToBlock( *sideIt );
 423       nbBlockSides--;
 424
 425       // edges of neighbour sides of B_FRONT corresponding to front's edges
 426       EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
 427       EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
 428       for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i )
 429         ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i]));
 430       if ( ok )
 431         ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP ));
 432
 433       // count the found sides
 434       _DUMP_(endl);
 435       for (int i = 0; i < B_UNDEFINED; ++i )
 436       {
 437         _DUMP_("** Block side "<< SBoxSides[i]);
 438         if ( block._side[ i ] )
 439         {
 440           if ( ok && i != B_FRONT)
 441           {
 442             setSideBoundToBlock( *block._side[ i ]._side );
 443             nbBlockSides--;
 444           }
 445           _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0));
 446           _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0));
 447           _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1));
 448           _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1));
 449         }
 450         else
 451         {
 452           _DUMP_("Not found"<<endl);
 453         }
 454       }
 455       if ( !ok )
 456         block._side[0] = block._side[1] = block._side[2] =
 457           block._side[3] = block._side[4] = block._side[5] = 0;
 458       else
 459         nbBlocks++;
 460     }
 461     if ( nbBlocks == 0 && _error.empty() )
 462       return BAD_MESH_ERR;
 463
 464     return nbBlocks;
 465   }
 466
 467   //================================================================================
 468   /*!
 469    * \brief Fill block side data starting from its corner quadrangle
 470    */
 471   //================================================================================
 472
 473   bool _Skin::fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad)
 474   {
 475     // Find out size of block side mesured in nodes and by the way find two rows
 476     // of nodes in two directions.
 477
 478     int x, y, nbX, nbY;
 479     const SMDS_MeshElement* firstQuad = cornerQuad;
 480     {
 481       // get corner node of cornerQuad
 482       const SMDS_MeshNode* nCorner = 0;
 483       for ( int i = 0; i < 4 && !nCorner; ++i )
 484         if ( isCornerNode( firstQuad->GetNode(i)))
 485           nCorner = firstQuad->GetNode(i);
 486       if ( !nCorner ) return false;
 487
 488       // get a node on block edge
 489       int iCorner = firstQuad->GetNodeIndex( nCorner );
 490       const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
 491
 492       // find out size of block side
 493       vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
 494       if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
 495            !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
 496         return false;
 497       nbX = horRow1.size(), nbY = verRow1.size();
 498
 499       // store found nodes
 500       side._index._xSize = horRow1.size();
 501       side._index._ySize = verRow1.size();
 502       side._grid.resize( side._index.size(), NULL );
 503
 504       for ( x = 0; x < horRow1.size(); ++x )
 505       {
 506         side.setNode( x, 0, horRow1[x] );
 507         side.setNode( x, 1, horRow2[x] );
 508       }
 509       for ( y = 0; y < verRow1.size(); ++y )
 510       {
 511         side.setNode( 0, y, verRow1[y] );
 512         side.setNode( 1, y, verRow2[y] );
 513       }
 514     }
 515     // Find the rest nodes
 516
 517     y = 1; // y of the row to fill
 518     TIDSortedElemSet emptySet, avoidSet;
 519     while ( ++y < nbY )
 520     {
 521       // get next firstQuad in the next row of quadrangles
 522       //
 523       //          n2up
 524       //     o---o               <- y row
 525       //     |   |
 526       //     o---o  o  o  o  o   <- found nodes
 527       //n1down    n2down
 528       //
 529       int i1down, i2down, i2up;
 530       const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
 531       const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
 532       avoidSet.clear(); avoidSet.insert( firstQuad );
 533       firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
 534                                                    &i1down, &i2down);
 535       if ( !isQuadrangle( firstQuad ))
 536         return BAD_MESH_ERR;
 537
 538       const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
 539       avoidSet.clear(); avoidSet.insert( firstQuad );
 540
 541       // find the rest nodes in the y-th row by faces in the row
 542
 543       x = 1;
 544       while ( ++x < nbX )
 545       {
 546         const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
 547                                                                         avoidSet, &i2up, &i2down);
 548         if ( !isQuadrangle( quad ))
 549           return BAD_MESH_ERR;
 550
 551         n2up   = oppositeNode( quad, i2down );
 552         n2down = oppositeNode( quad, i2up );
 553         avoidSet.clear(); avoidSet.insert( quad );
 554
 555         side.setNode( x, y, n2up );
 556       }
 557     }
 558
 559     return true;
 560   }
 561
 562   //================================================================================
 563   /*!
 564    * \brief Try to find a block side adjacent to the given side by given edge
 565    */
 566   //================================================================================
 567
 568   _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
 569                                            EQuadEdge sharedSideEdge1,
 570                                            EQuadEdge sharedSideEdge2)
 571   {
 572     _Block& block = _blocks.back();
 573     _OrientedBlockSide& side1 = block._side[ startBlockSide ];
 574
 575     // get corner nodes of the given block edge
 576     bool xMax1=0, yMax1=0, xMax2=1, yMax2=1;
 577     switch( sharedSideEdge1 )
 578     {
 579     case Q_BOTTOM: yMax2 = 0; break;
 580     case Q_RIGHT:  xMax1 = 1; break;
 581     case Q_TOP:    yMax1 = 1; break;
 582     case Q_LEFT:   xMax2 = 0; break;
 583     default:
 584       error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__);
 585       return 0;
 586     }
 587     const SMDS_MeshNode* n1 = side1.cornerNode( xMax1, yMax1);
 588     const SMDS_MeshNode* n2 = side1.cornerNode( xMax2, yMax2);
 589
 590     set< _BlockSide* >& sidesWithN1 = _corner2sides[ n1 ];
 591     set< _BlockSide* >& sidesWithN2 = _corner2sides[ n2 ];
 592     if ( sidesWithN1.empty() || sidesWithN2.empty() )
 593     {
 594       _DUMP_("no sides by nodes "<< n1->GetID() << " " << n2->GetID() );
 595       return 0;
 596     }
 597
 598     // find all sides sharing both nodes n1 and n2
 599     set< _BlockSide* > sidesOnEdge;
 600     set< _BlockSide* >::iterator sideIt;
 601     std::set_intersection( sidesWithN1.begin(), sidesWithN1.end(),
 602                            sidesWithN2.begin(), sidesWithN2.end(),
 603                            inserter( sidesOnEdge, sidesOnEdge.begin()));
 604
 605     // exclude loaded sides of block from sidesOnEdge
 606     int nbLoadedSides = 0;
 607     for (int i = 0; i < B_UNDEFINED; ++i )
 608     {
 609       if ( block._side[ i ] )
 610       {
 611         nbLoadedSides++;
 612         sidesOnEdge.erase( block._side[ i ]._side );
 613       }
 614     }
 615     int nbSidesOnEdge = sidesOnEdge.size();
 616     _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge);
 617     if ( nbSidesOnEdge == 0 )
 618       return 0;
 619
 620     _BlockSide* foundSide = 0;
 621     if ( nbSidesOnEdge == 1 || nbSidesOnEdge == 2 && nbLoadedSides == 1 )
 622     {
 623       foundSide = *sidesOnEdge.begin();
 624     }
 625     else
 626     {
 627       // Select one of found sides most close to startBlockSide
 628
 629       // gravity center of already loaded block sides
 630       gp_XYZ gc;
 631       for (int i = 0; i < B_UNDEFINED; ++i )
 632         if ( block._side[ i ] )
 633           gc += block._side[ i ]._side->getGC();
 634       gc /= nbLoadedSides;
 635
 636       gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
 637       gp_Vec p2p1( p1 - p2 );
 638
 639       const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
 640       gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
 641       gp_Vec side1Dir( p1, p1Op );
 642
 643       map < double , _BlockSide* > angleOfSide;
 644       for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
 645       {
 646         _BlockSide* sideI = *sideIt;
 647         const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
 648         gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
 649         gp_Vec sideIDir( p1, p1Op );
 650         double angle = side1Dir.AngleWithRef( sideIDir, p2p1 );
 651         if ( angle < 0 ) angle += 2 * PI;
 652         angleOfSide.insert( make_pair( angle, sideI ));
 653       }
 654       gp_Vec gcDir( p1, gc );
 655       double gcAngle = side1Dir.AngleWithRef( gcDir, p2p1 );
 656       foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
 657     }
 658
 659     // Orient the found side correctly
 660
 661     // corners of found side corresponding to nodes n1 and n2
 662     xMax1= yMax1 = 0; xMax2 = yMax2 = 1;
 663     switch( sharedSideEdge2 )
 664     {
 665     case Q_BOTTOM: yMax2 = 0; break;
 666     case Q_RIGHT:  xMax1 = 1; break;
 667     case Q_TOP:    yMax1 = 1; break;
 668     case Q_LEFT:   xMax2 = 0; break;
 669     default:
 670       error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__);
 671       return 0;
 672     }
 673     for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
 674     {
 675       _OrientedBlockSide orientedSide( foundSide, ori );
 676       const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
 677       const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
 678       if ( n1 == n12 && n2 == n22 )
 679         return orientedSide;
 680     }
 681     error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
 682           << " of side " << startBlockSide
 683           << " of block " << _blocks.size());
 684     return 0;
 685   }
 686
 687   //================================================================================
 688   /*!
 689    * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
 690    * from the given quadrangle until another block corner encounters.
 691    *  n1 and n2 are at bottom of quad, n1 is at block corner.
 692    */
 693   //================================================================================
 694
 695   bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement*       quad,
 696                                   const SMDS_MeshNode*          n1,
 697                                   const SMDS_MeshNode*          n2,
 698                                   vector<const SMDS_MeshNode*>& row1,
 699                                   vector<const SMDS_MeshNode*>& row2,
 700                                   const bool                    alongN1N2 )
 701   {
 702     const SMDS_MeshNode* corner1 = n1;
 703
 704     // Store nodes of quad in the rows and find new n1 and n2 to get
 705     // the next face so that new n2 is on block edge
 706     int i1 = quad->GetNodeIndex( n1 );
 707     int i2 = quad->GetNodeIndex( n2 );
 708     row1.clear(); row2.clear();
 709     row1.push_back( n1 );
 710     if ( alongN1N2 )
 711     {
 712       row1.push_back( n2 );
 713       row2.push_back( oppositeNode( quad, i2 ));
 714       row2.push_back( n1 = oppositeNode( quad, i1 ));
 715     }
 716     else
 717     {
 718       row2.push_back( n2 );
 719       row1.push_back( n2 = oppositeNode( quad, i2 ));
 720       row2.push_back( n1 = oppositeNode( quad, i1 ));
 721     }
 722
 723     // Find the rest nodes
 724     TIDSortedElemSet emptySet, avoidSet;
 725     while ( !isCornerNode( n2 ))
 726     {
 727       avoidSet.clear(); avoidSet.insert( quad );
 728       quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
 729       if ( !isQuadrangle( quad ))
 730         return BAD_MESH_ERR;
 731
 732       row1.push_back( n2 = oppositeNode( quad, i1 ));
 733       row2.push_back( n1 = oppositeNode( quad, i2 ));
 734     }
 735     return n1 != corner1;
 736   }
 737
 738   //================================================================================
 739   /*!
 740    * \brief Return a corner face by a corner node
 741    */
 742   //================================================================================
 743
 744   const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
 745   {
 746     int x, y, isXMax, isYMax, found = 0;
 747     for ( isXMax = 0; isXMax < 2; ++isXMax )
 748     {
 749       for ( isYMax = 0; isYMax < 2; ++isYMax )
 750       {
 751         x = isXMax ? _index._xSize-1 : 0;
 752         y = isYMax ? _index._ySize-1 : 0;
 753         found = ( getNode(x,y) == cornerNode );
 754         if ( found ) break;
 755       }
 756       if ( found ) break;
 757     }
 758     if ( !found ) return 0;
 759     int dx = isXMax ? -1 : +1;
 760     int dy = isYMax ? -1 : +1;
 761     const SMDS_MeshNode* n1 = getNode(x,y);
 762     const SMDS_MeshNode* n2 = getNode(x+dx,y);
 763     const SMDS_MeshNode* n3 = getNode(x,y+dy);
 764     const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
 765     return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
 766   }
 767
 768 }
 769
 770 //=======================================================================
 771 //function : StdMeshers_HexaFromSkin_3D
 772 //purpose  :
 773 //=======================================================================
 774
 775 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen)
 776   :SMESH_3D_Algo(hypId, studyId, gen)
 777 {
 778   MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D");
 779   _name = "HexaFromSkin_3D";
 780 }
 781
 782 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
 783 {
 784   MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D");
 785 }
 786
 787 //================================================================================
 788 /*!
 789  * \brief Main method, which generates hexaheda
 790  */
 791 //================================================================================
 792
 793 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
 794 {
 795   _Skin skin;
 796   int nbBlocks = skin.findBlocks(aMesh);
 797   if ( nbBlocks == 0 )
 798     return error( skin.error());
 799
 800   vector< vector< const SMDS_MeshNode* > > columns;
 801   int x, xSize, y, ySize, z, zSize;
 802   _Indexer colIndex;
 803
 804   for ( int i = 0; i < nbBlocks; ++i )
 805   {
 806     const _Block& block = skin.getBlock( i );
 807
 808     // ------------------------------------------
 809     // Fill columns of nodes with existing nodes
 810     // ------------------------------------------
 811
 812     xSize = block.getSide(B_BOTTOM).getHoriSize();
 813     ySize = block.getSide(B_BOTTOM).getVertSize();
 814     zSize = block.getSide(B_FRONT ).getVertSize();
 815     int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
 816     colIndex = _Indexer( xSize, ySize );
 817     columns.resize( colIndex.size() );
 818
 819     // fill node columns by front and back box sides
 820     for ( x = 0; x < xSize; ++x ) {
 821       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
 822       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
 823       column0.resize( zSize );
 824       column1.resize( zSize );
 825       for ( z = 0; z < zSize; ++z ) {
 826         column0[ z ] = block.getSide(B_FRONT).node( x, z );
 827         column1[ z ] = block.getSide(B_BACK) .node( x, z );
 828       }
 829     }
 830     // fill node columns by left and right box sides
 831     for ( y = 1; y < ySize-1; ++y ) {
 832       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
 833       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
 834       column0.resize( zSize );
 835       column1.resize( zSize );
 836       for ( z = 0; z < zSize; ++z ) {
 837         column0[ z ] = block.getSide(B_LEFT) .node( y, z );
 838         column1[ z ] = block.getSide(B_RIGHT).node( y, z );
 839       }
 840     }
 841     // get nodes from top and bottom box sides
 842     for ( x = 1; x < xSize-1; ++x ) {
 843       for ( y = 1; y < ySize-1; ++y ) {
 844         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 845         column.resize( zSize );
 846         column.front() = block.getSide(B_BOTTOM).node( x, y );
 847         column.back()  = block.getSide(B_TOP)   .node( x, y );
 848       }
 849     }
 850
 851     // ----------------------------
 852     // Add internal nodes of a box
 853     // ----------------------------
 854     // projection points of internal nodes on box subshapes by which
 855     // coordinates of internal nodes are computed
 856     vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
 857
 858     // projections on vertices are constant
 859     pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
 860     pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
 861     pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
 862     pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
 863     pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
 864     pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
 865     pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
 866     pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
 867
 868     for ( x = 1; x < xSize-1; ++x )
 869     {
 870       gp_XYZ params; // normalized parameters of internal node within a unit box
 871       params.SetCoord( 1, x / double(X) );
 872       for ( y = 1; y < ySize-1; ++y )
 873       {
 874         params.SetCoord( 2, y / double(Y) );
 875         // column to fill during z loop
 876         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 877         // projections on horizontal edges
 878         pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
 879         pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
 880         pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
 881         pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
 882         pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
 883         pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
 884         pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
 885         pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
 886         // projections on horizontal sides
 887         pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
 888         pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP)   .xyz( x, y );
 889         for ( z = 1; z < zSize-1; ++z ) // z loop
 890         {
 891           params.SetCoord( 3, z / double(Z) );
 892           // projections on vertical edges
 893           pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
 894           pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
 895           pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
 896           pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
 897           // projections on vertical sides
 898           pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
 899           pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
 900           pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
 901           pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
 902
 903           // compute internal node coordinates
 904           gp_XYZ coords;
 905           SMESH_Block::ShellPoint( params, pointOnShape, coords );
 906           column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
 907
 908 #ifdef DEB_GRID
 909           // debug
 910           //cout << "----------------------------------------------------------------------"<<endl;
 911           //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
 912           //{
 913           //  gp_XYZ p = pointOnShape[ id ];
 914           //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
 915           //}
 916           //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
 917           //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
 918 #endif
 919         }
 920       }
 921     }
 922     // ----------------
 923     // Add hexahedrons
 924     // ----------------
 925
 926     // find out orientation
 927     const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
 928     const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
 929     const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
 930     const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
 931     const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
 932     const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
 933     const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
 934     const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
 935     SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
 936                                     n001,n011,n111,n101);
 937     bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
 938
 939     // add elements
 940     for ( x = 0; x < xSize-1; ++x ) {
 941       for ( y = 0; y < ySize-1; ++y ) {
 942         vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
 943         vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
 944         vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
 945         vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
 946         // bottom face normal of a hexa mush point outside the volume
 947         if ( isForw )
 948           for ( z = 0; z < zSize-1; ++z )
 949             aHelper->AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
 950                                col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
 951         else
 952           for ( z = 0; z < zSize-1; ++z )
 953             aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
 954                                col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
 955       }
 956     }
 957   } // loop on blocks
 958
 959   return true;
 960 }
 961
 962 //================================================================================
 963 /*!
 964  * \brief Evaluate nb of hexa
 965  */
 966 //================================================================================
 967
 968 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh &         aMesh,
 969                                           const TopoDS_Shape & aShape,
 970                                           MapShapeNbElems&     aResMap)
 971 {
 972   _Skin skin;
 973   int nbBlocks = skin.findBlocks(aMesh);
 974   if ( nbBlocks == 0 )
 975     return error( skin.error());
 976
 977   bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
 978
 979   int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
 980   vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
 981   if ( entity >= nbByType.size() )
 982     nbByType.resize( SMDSEntity_Last, 0 );
 983
 984   for ( int i = 0; i < nbBlocks; ++i )
 985   {
 986     const _Block& block = skin.getBlock( i );
 987
 988     int nbX = block.getSide(B_BOTTOM).getHoriSize();
 989     int nbY = block.getSide(B_BOTTOM).getVertSize();
 990     int nbZ = block.getSide(B_FRONT ).getVertSize();
 991
 992     int nbHexa  = (nbX-1) * (nbY-1) * (nbZ-1);
 993     int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
 994     if ( secondOrder )
 995       nbNodes +=
 996         (nbX-2) * (nbY-2) * (nbZ-1) +
 997         (nbX-2) * (nbY-1) * (nbZ-2) +
 998         (nbX-1) * (nbY-2) * (nbZ-2);
 999
1000
1001     nbByType[ entity ] += nbHexa;
1002     nbByType[ SMDSEntity_Node ] += nbNodes;
1003   }
1004
1005   return true;
1006 }
1007
1008 //================================================================================
1009 /*!
1010  * \brief Abstract method must be defined but does nothing
1011  */
1012 //================================================================================
1013
1014 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1015                                                  Hypothesis_Status& aStatus)
1016 {
1017   aStatus = SMESH_Hypothesis::HYP_OK;
1018   return true;
1019 }
1020
1021 //================================================================================
1022 /*!
1023  * \brief Abstract method must be defined but just reports an error as this
1024  *  algo is not intended to work with shapes
1025  */
1026 //================================================================================
1027
1028 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1029 {
1030   return error("Algorithm can't work with geometrical shapes");
1031 }
1032