src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx

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