src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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