src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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