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