src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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