src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx

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