src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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