src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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