src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File      : StdMeshers_HexaFromSkin_3D.cxx
  21 // Created   : Wed Jan 27 12:28:07 2010
  22 // Author    : Edward AGAPOV (eap)
  23
  24 #include "StdMeshers_HexaFromSkin_3D.hxx"
  25
  26 #include "SMDS_VolumeOfNodes.hxx"
  27 #include "SMDS_VolumeTool.hxx"
  28 #include "SMESHDS_Mesh.hxx"
  29 #include "SMESH_Block.hxx"
  30 #include "SMESH_Mesh.hxx"
  31 #include "SMESH_MeshAlgos.hxx"
  32 #include "SMESH_MesherHelper.hxx"
  33
  34 #include <gp_Ax2.hxx>
  35
  36 #include <limits>
  37
  38 using namespace std;
  39
  40 // Define error message and _MYDEBUG_ if needed
  41 #ifdef _DEBUG_
  42 #define BAD_MESH_ERR \
  43   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
  44                       __FILE__ ":" )<<__LINE__)
  45 //#define _MYDEBUG_
  46 #else
  47 #define BAD_MESH_ERR \
  48   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
  49 #endif
  50
  51
  52 // Debug output
  53 #ifdef _MYDEBUG_
  54 #define _DUMP_(msg) cout << msg << endl
  55 #else
  56 #define _DUMP_(msg)
  57 #endif
  58
  59
  60 namespace
  61 {
  62   enum EBoxSides //!< sides of the block
  63     {
  64       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
  65     };
  66 #ifdef _MYDEBUG_
  67   const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only
  68     {
  69       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
  70     };
  71 #endif
  72   enum EQuadEdge //!< edges of quadrangle side
  73     {
  74       Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
  75     };
  76
  77
  78   //================================================================================
  79   /*!
  80    * \brief return logical coordinates (i.e. min or max) of ends of edge
  81    */
  82   //================================================================================
  83
  84   bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
  85   {
  86     xMax1=0, yMax1=0, xMax2=1, yMax2=1;
  87     switch( edge )
  88     {
  89     case Q_BOTTOM: yMax2 = 0; break;
  90     case Q_RIGHT:  xMax1 = 1; break;
  91     case Q_TOP:    yMax1 = 1; break;
  92     case Q_LEFT:   xMax2 = 0; break;
  93     default:
  94       return false;
  95     }
  96     return true;
  97   }
  98
  99   //================================================================================
 100   /*!
 101    * \brief return true if a node is at block corner
 102    *
 103    * This check is valid for simple cases only
 104    */
 105   //================================================================================
 106
 107   bool isCornerNode( const SMDS_MeshNode* n )
 108   {
 109     int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
 110     if ( nbF % 2 )
 111       return true;
 112
 113     set<const SMDS_MeshNode*> nodesInInverseFaces;
 114     SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
 115     while ( fIt->more() )
 116     {
 117       const SMDS_MeshElement* face = fIt->next();
 118       nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
 119     }
 120
 121     return (int)nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
 122   }
 123
 124   //================================================================================
 125   /*!
 126    * \brief check element type
 127    */
 128   //================================================================================
 129
 130   bool isQuadrangle(const SMDS_MeshElement* e)
 131   {
 132     return ( e && e->NbCornerNodes() == 4 );
 133   }
 134
 135   //================================================================================
 136   /*!
 137    * \brief return opposite node of a quadrangle face
 138    */
 139   //================================================================================
 140
 141   const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
 142   {
 143     return quad->GetNode( (iNode+2) % 4 );
 144   }
 145
 146   //================================================================================
 147   /*!
 148    * \brief Converter of a pair of integers to a sole index
 149    */
 150   struct _Indexer
 151   {
 152     int _xSize, _ySize;
 153     _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
 154     int size() const { return _xSize * _ySize; }
 155     int operator()(int x, int y) const { return y * _xSize + x; }
 156   };
 157   //================================================================================
 158   /*!
 159    * \brief Oriented converter of a pair of integers to a sole index
 160    */
 161   class _OrientedIndexer : public _Indexer
 162   {
 163   public:
 164     enum OriFlags //!< types of block side orientation
 165       {
 166         REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
 167       };
 168     _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
 169       _Indexer( indexer._xSize, indexer._ySize ),
 170       _xSize (indexer._xSize), _ySize(indexer._ySize),
 171       _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
 172       _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
 173       _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
 174     {
 175       (*_swapFun)( _xSize, _ySize );
 176     }
 177     //!< Return index by XY
 178     int operator()(int x, int y) const
 179     {
 180       (*_xRevFun)( x, const_cast<int&>( _xSize ));
 181       (*_yRevFun)( y, const_cast<int&>( _ySize ));
 182       (*_swapFun)( x, y );
 183       return _Indexer::operator()(x,y);
 184     }
 185     //!< Return index for a corner
 186     int corner(bool xMax, bool yMax) const
 187     {
 188       int x = xMax, y = yMax, size = 2;
 189       (*_xRevFun)( x, size );
 190       (*_yRevFun)( y, size );
 191       (*_swapFun)( x, y );
 192       return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
 193     }
 194     int xSize() const { return _xSize; }
 195     int ySize() const { return _ySize; }
 196   private:
 197     _Indexer _indexer;
 198     int _xSize, _ySize;
 199
 200     typedef void (*TFun)(int& x, int& y);
 201     TFun _xRevFun, _yRevFun, _swapFun;
 202
 203     static void lazy   (int&, int&) {}
 204     static void reverse(int& x, int& size) { x = size - x - 1; }
 205     static void swap   (int& x, int& y) { std::swap(x,y); }
 206   };
 207   //================================================================================
 208   /*!
 209    * \brief Structure corresponding to the meshed side of block
 210    */
 211   struct _BlockSide
 212   {
 213     vector<const SMDS_MeshNode*> _grid;
 214     _Indexer                     _index;
 215     int                          _nbBlocksExpected;
 216     int                          _nbBlocksFound;
 217
 218 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
 219 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < (int)pobj->_grid.size()) ? i : int(1e100)]
 220 #else
 221 #define _grid_access_(pobj, i) pobj->_grid[ i ]
 222 #endif
 223     //!< Return node at XY
 224     const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
 225     //!< Set node at XY
 226     void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
 227     //!< Return an edge
 228     SMESH_OrientedLink getEdge(EQuadEdge edge) const
 229     {
 230       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
 231       return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
 232     }
 233     //!< Return a corner node
 234     const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
 235     {
 236       return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
 237     }
 238     const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
 239     //!< True if all blocks this side belongs to have been found
 240     bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
 241     //!< Return coordinates of node at XY
 242     gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
 243     //!< Return gravity center of the four corners and the middle node
 244     gp_XYZ getGC() const
 245     {
 246       gp_XYZ xyz =
 247         getXYZ( 0, 0 ) +
 248         getXYZ( _index._xSize-1, 0 ) +
 249         getXYZ( 0, _index._ySize-1 ) +
 250         getXYZ( _index._xSize-1, _index._ySize-1 ) +
 251         getXYZ( _index._xSize/2, _index._ySize/2 );
 252       return xyz / 5;
 253     }
 254     //!< Return number of mesh faces
 255     int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
 256   };
 257   //================================================================================
 258   /*!
 259    * \brief _BlockSide with changed orientation
 260    */
 261   struct _OrientedBlockSide
 262   {
 263     _BlockSide*       _side;
 264     _OrientedIndexer  _index;
 265
 266     _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
 267       _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
 268     //!< return coordinates by XY
 269     gp_XYZ xyz(int x, int y) const
 270     {
 271       return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
 272     }
 273     //!< safely return a node by XY
 274     const SMDS_MeshNode* node(int x, int y) const
 275     {
 276       size_t i = _index( x, y );
 277       return ( i >= _side->_grid.size() ) ? 0 : _side->_grid[i];
 278     }
 279     //!< Return an edge
 280     SMESH_OrientedLink edge(EQuadEdge edge) const
 281     {
 282       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
 283       return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
 284     }
 285     //!< Return a corner node
 286     const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
 287     {
 288       return _grid_access_(_side, _index.corner( isXMax, isYMax ));
 289     }
 290     //!< return its size in nodes
 291     int getHoriSize() const { return _index.xSize(); }
 292     int getVertSize() const  { return _index.ySize(); }
 293     //!< True if _side has been initialized
 294     operator bool() const { return _side; }
 295     //! Direct access to _side
 296     const _BlockSide* operator->() const { return _side; }
 297     _BlockSide* operator->() { return _side; }
 298   };
 299   //================================================================================
 300   /*!
 301    * \brief Meshed skin of block
 302    */
 303   struct _Block
 304   {
 305     _OrientedBlockSide        _side[6]; // 6 sides of a sub-block
 306     set<const SMDS_MeshNode*> _corners;
 307
 308     const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
 309     bool setSide( int i, const _OrientedBlockSide& s)
 310     {
 311       if (( _side[i] = s ))
 312       {
 313         _corners.insert( s.cornerNode(0,0));
 314         _corners.insert( s.cornerNode(1,0));
 315         _corners.insert( s.cornerNode(0,1));
 316         _corners.insert( s.cornerNode(1,1));
 317       }
 318       return s;
 319     }
 320     void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
 321     bool hasSide( const _OrientedBlockSide& s) const
 322     {
 323       if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
 324       return false;
 325     }
 326     int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
 327     bool isValid() const;
 328   };
 329   //================================================================================
 330   /*!
 331    * \brief Skin mesh possibly containing several meshed blocks
 332    */
 333   class _Skin
 334   {
 335   public:
 336
 337     int findBlocks(SMESH_Mesh& mesh);
 338     //!< return i-th block
 339     const _Block& getBlock(int i) const { return _blocks[i]; }
 340     //!< return error description
 341     const SMESH_Comment& error() const { return _error; }
 342
 343   private:
 344     bool fillSide( _BlockSide&             side,
 345                    const SMDS_MeshElement* cornerQuad,
 346                    const SMDS_MeshNode*    cornerNode);
 347     bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
 348                              const SMDS_MeshNode*    n1,
 349                              const SMDS_MeshNode*    n2,
 350                              vector<const SMDS_MeshNode*>& verRow1,
 351                              vector<const SMDS_MeshNode*>& verRow2,
 352                              bool alongN1N2 );
 353     _OrientedBlockSide findBlockSide( EBoxSides           startBlockSide,
 354                                       EQuadEdge           sharedSideEdge1,
 355                                       EQuadEdge           sharedSideEdge2,
 356                                       bool                withGeometricAnalysis,
 357                                       set< _BlockSide* >& sidesAround);
 358     //!< update own data and data of the side bound to block
 359     void setSideBoundToBlock( _BlockSide& side )
 360     {
 361       if ( side._nbBlocksFound++, side.isBound() )
 362         for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
 363           _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
 364     }
 365     //!< store reason of error
 366     int error(const SMESH_Comment& reason) { _error = reason; return 0; }
 367
 368     SMESH_Comment      _error;
 369
 370     list< _BlockSide > _allSides;
 371     vector< _Block >   _blocks;
 372
 373     //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
 374     map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
 375   };
 376
 377   //================================================================================
 378   /*!
 379    * \brief Find blocks and return their number
 380    */
 381   //================================================================================
 382
 383   int _Skin::findBlocks(SMESH_Mesh& mesh)
 384   {
 385     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 386
 387     // Find a node at any block corner
 388
 389     SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
 390     if ( !nIt->more() ) return error("Empty mesh");
 391
 392     const SMDS_MeshNode* nCorner = 0;
 393     while ( nIt->more() )
 394     {
 395       nCorner = nIt->next();
 396       if ( isCornerNode( nCorner ))
 397         break;
 398       else
 399         nCorner = 0;
 400     }
 401     if ( !nCorner )
 402       return BAD_MESH_ERR;
 403
 404     // --------------------------------------------------------------------
 405     // Find all block sides starting from mesh faces sharing the corner node
 406     // --------------------------------------------------------------------
 407
 408     int nbFacesOnSides = 0;
 409     TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
 410     list< const SMDS_MeshNode* > corners( 1, nCorner );
 411     list< const SMDS_MeshNode* >::iterator corner = corners.begin();
 412     while ( corner != corners.end() )
 413     {
 414       SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
 415       while ( faceIt->more() )
 416       {
 417         const SMDS_MeshElement* face = faceIt->next();
 418         if ( !cornerFaces.insert( face ).second )
 419           continue; // already loaded block side
 420
 421         if ( !isQuadrangle( face ))
 422           return error("Non-quadrangle elements in the input mesh");
 423
 424         if ( _allSides.empty() || !_allSides.back()._grid.empty() )
 425           _allSides.push_back( _BlockSide() );
 426
 427         _BlockSide& side = _allSides.back();
 428         if ( !fillSide( side, face, *corner ))
 429         {
 430           if ( !_error.empty() )
 431             return false;
 432         }
 433         else
 434         {
 435           for ( int isXMax = 0; isXMax < 2; ++isXMax )
 436             for ( int isYMax = 0; isYMax < 2; ++isYMax )
 437             {
 438               const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
 439               corners.push_back( nCorner );
 440               cornerFaces.insert( side.getCornerFace( nCorner ));
 441             }
 442           for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
 443             _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
 444
 445           nbFacesOnSides += side.getNbFaces();
 446         }
 447       }
 448       ++corner;
 449
 450       // find block sides of other domains if any
 451       if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
 452       {
 453         while ( nIt->more() )
 454         {
 455           nCorner = nIt->next();
 456           if ( isCornerNode( nCorner ))
 457             corner = corners.insert( corner, nCorner );
 458         }
 459         nbFacesOnSides = mesh.NbQuadrangles();
 460       }
 461     }
 462
 463     if ( _allSides.empty() )
 464       return BAD_MESH_ERR;
 465     if ( _allSides.back()._grid.empty() )
 466       _allSides.pop_back();
 467     _DUMP_("Nb detected sides "<< _allSides.size());
 468
 469     // ---------------------------
 470     // Organize sides into blocks
 471     // ---------------------------
 472
 473     // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
 474     int nbBlockSides = 0; // total nb of block sides taking into account their sharing
 475     multimap<int, _BlockSide* > sortedSides;
 476     {
 477       list < _BlockSide >::iterator sideIt = _allSides.begin();
 478       for ( ; sideIt != _allSides.end(); ++sideIt )
 479       {
 480         _BlockSide& side = *sideIt;
 481         bool isSharedSide = true;
 482         int nbAdjacent = 0;
 483         for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
 484         {
 485           int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
 486           nbAdjacent += nbAdj;
 487           isSharedSide = ( nbAdj > 2 );
 488         }
 489         side._nbBlocksFound = 0;
 490         side._nbBlocksExpected = isSharedSide ? 2 : 1;
 491         nbBlockSides += side._nbBlocksExpected;
 492         sortedSides.insert( make_pair( nbAdjacent, & side ));
 493       }
 494     }
 495
 496     // find sides of each block
 497     int nbBlocks = 0;
 498     while ( nbBlockSides >= 6 )
 499     {
 500       // get any side not bound to all blocks it belongs to
 501       multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
 502       while ( i_side != sortedSides.end() && i_side->second->isBound())
 503         ++i_side;
 504
 505       // start searching for block sides from the got side
 506       bool ok = true;
 507       if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
 508         _blocks.resize( _blocks.size() + 1 );
 509
 510       _Block& block = _blocks.back();
 511       block.setSide( B_FRONT, i_side->second );
 512       setSideBoundToBlock( *i_side->second );
 513       nbBlockSides--;
 514
 515       // edges of adjacent sides of B_FRONT corresponding to front's edges
 516       EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
 517       EQuadEdge edgeOfAdj  [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
 518       // first find all sides detectable w/o advanced analysis,
 519       // then repeat the search, which then may pass without advanced analysis
 520       set< _BlockSide* > sidesAround;
 521       for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
 522       {
 523         // try to find 4 sides adjacent to a FRONT side
 524         for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
 525           if ( !block._side[i] )
 526             ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
 527                                                   advAnalys, sidesAround));
 528         // try to find a BACK side by a TOP one
 529         if ( ok || !advAnalys )
 530           if ( !block._side[B_BACK] && block._side[B_TOP] )
 531             ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
 532                                                        advAnalys, sidesAround ));
 533         if ( !advAnalys ) ok = true;
 534       }
 535       ok = block.isValid();
 536       if ( ok )
 537       {
 538         // check if just found block is same as one of previously found blocks
 539         bool isSame = false;
 540         for ( size_t i = 1; i < _blocks.size() && !isSame; ++i )
 541           isSame = ( block._corners == _blocks[i-1]._corners );
 542         ok = !isSame;
 543       }
 544
 545       // count the found sides
 546       _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
 547       for (int i = 0; i < NB_BLOCK_SIDES; ++i )
 548       {
 549         _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
 550         if ( block._side[ i ] )
 551         {
 552           if ( ok && i != B_FRONT)
 553           {
 554             setSideBoundToBlock( *block._side[ i ]._side );
 555             nbBlockSides--;
 556           }
 557           _DUMP_("\t corners "<<
 558                  block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
 559                  block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
 560                  block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
 561                  block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
 562         }
 563         else
 564         {
 565           _DUMP_("\t not found"<<endl);
 566         }
 567       }
 568       if ( !ok )
 569         block.clear();
 570       else
 571         nbBlocks++;
 572     }
 573     _DUMP_("Nb found blocks "<< nbBlocks <<endl);
 574
 575     if ( nbBlocks == 0 && _error.empty() )
 576       return BAD_MESH_ERR;
 577
 578     return nbBlocks;
 579   }
 580
 581   //================================================================================
 582   /*!
 583    * \brief Fill block side data starting from its corner quadrangle
 584    */
 585   //================================================================================
 586
 587   bool _Skin::fillSide( _BlockSide&             side,
 588                         const SMDS_MeshElement* cornerQuad,
 589                         const SMDS_MeshNode*    nCorner)
 590   {
 591     // Find out size of block side measured in nodes and by the way find two rows
 592     // of nodes in two directions.
 593
 594     int x, y, nbX, nbY;
 595     const SMDS_MeshElement* firstQuad = cornerQuad;
 596     {
 597       // get a node on block edge
 598       int iCorner = firstQuad->GetNodeIndex( nCorner );
 599       const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
 600
 601       // find out size of block side
 602       vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
 603       if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
 604            !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
 605         return false;
 606       nbX = horRow1.size(), nbY = verRow1.size();
 607
 608       // store found nodes
 609       side._index._xSize = horRow1.size();
 610       side._index._ySize = verRow1.size();
 611       side._grid.resize( side._index.size(), NULL );
 612
 613       for ( x = 0; x < nbX; ++x )
 614       {
 615         side.setNode( x, 0, horRow1[x] );
 616         side.setNode( x, 1, horRow2[x] );
 617       }
 618       for ( y = 0; y < nbY; ++y )
 619       {
 620         side.setNode( 0, y, verRow1[y] );
 621         side.setNode( 1, y, verRow2[y] );
 622       }
 623     }
 624     // Find the rest nodes
 625
 626     y = 1; // y of the row to fill
 627     TIDSortedElemSet emptySet, avoidSet;
 628     while ( ++y < nbY )
 629     {
 630       // get next firstQuad in the next row of quadrangles
 631       //
 632       //          n2up
 633       //     o---o               <- y row
 634       //     |   |
 635       //     o---o  o  o  o  o   <- found nodes
 636       //n1down    n2down
 637       //
 638       int i1down, i2down, i2up;
 639       const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
 640       const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
 641       avoidSet.clear(); avoidSet.insert( firstQuad );
 642       firstQuad = SMESH_MeshAlgos::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
 643                                                    &i1down, &i2down);
 644       if ( !isQuadrangle( firstQuad ))
 645         return BAD_MESH_ERR;
 646
 647       const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
 648       avoidSet.clear(); avoidSet.insert( firstQuad );
 649
 650       // find the rest nodes in the y-th row by faces in the row
 651
 652       x = 1;
 653       while ( ++x < nbX )
 654       {
 655         const SMDS_MeshElement* quad = SMESH_MeshAlgos::FindFaceInSet( n2up, n2down, emptySet,
 656                                                                         avoidSet, &i2up, &i2down);
 657         if ( !isQuadrangle( quad ))
 658           return BAD_MESH_ERR;
 659
 660         n2up   = oppositeNode( quad, i2down );
 661         n2down = oppositeNode( quad, i2up );
 662         avoidSet.clear(); avoidSet.insert( quad );
 663
 664         side.setNode( x, y, n2up );
 665       }
 666     }
 667
 668     // check side validity
 669     bool ok =
 670       side.getCornerFace( side.getCornerNode( 0, 0 )) &&
 671       side.getCornerFace( side.getCornerNode( 1, 0 )) &&
 672       side.getCornerFace( side.getCornerNode( 0, 1 )) &&
 673       side.getCornerFace( side.getCornerNode( 1, 1 ));
 674
 675     return ok;
 676   }
 677
 678   //================================================================================
 679   /*!
 680    * \brief Return true if it's possible to make a loop over corner2Sides starting
 681    * from the startSide
 682    */
 683   //================================================================================
 684
 685   bool isClosedChainOfSides( _BlockSide*                                        startSide,
 686                              map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
 687   {
 688     // get start and end nodes
 689     const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
 690     for ( int y = 0; y < 2; ++y )
 691       for ( int x = 0; x < 2; ++x )
 692       {
 693         n = startSide->getCornerNode(x,y);
 694         if ( !corner2Sides.count( n )) continue;
 695         if ( n1 )
 696           n2 = n;
 697         else
 698           n1 = n;
 699       }
 700     if ( !n2 ) return false;
 701
 702     map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
 703       c2sides = corner2Sides.find( n1 );
 704     if ( c2sides == corner2Sides.end() ) return false;
 705
 706     int nbChainLinks = 1;
 707     n = n1;
 708     _BlockSide* prevSide = startSide;
 709     while ( n != n2 )
 710     {
 711       // get the next side sharing n
 712       list< _BlockSide* > & sides = c2sides->second;
 713       _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
 714       if ( nextSide == prevSide ) return false;
 715
 716       // find the next corner of the nextSide being in corner2Sides
 717       n1 = n;
 718       n = 0;
 719       for ( int y = 0; y < 2 && !n; ++y )
 720         for ( int x = 0; x < 2; ++x )
 721         {
 722           n = nextSide->getCornerNode(x,y);
 723           c2sides = corner2Sides.find( n );
 724           if ( n == n1 || c2sides == corner2Sides.end() )
 725             n = 0;
 726           else
 727             break;
 728         }
 729       if ( !n ) return false;
 730
 731       prevSide = nextSide;
 732
 733       if ( ++nbChainLinks > NB_QUAD_SIDES )
 734         return false;
 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_MeshAlgos::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     ok = ok && ( _side[ B_BOTTOM ]._index.size() == _side[ B_TOP  ]._index.size() &&
1013                  _side[ B_RIGHT  ]._index.size() == _side[ B_LEFT ]._index.size() &&
1014                  _side[ B_FRONT  ]._index.size() == _side[ B_BACK ]._index.size() );
1015     return ok;
1016   }
1017
1018 } // namespace
1019
1020 //=======================================================================
1021 //function : StdMeshers_HexaFromSkin_3D
1022 //purpose  :
1023 //=======================================================================
1024
1025 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, SMESH_Gen* gen)
1026   :SMESH_3D_Algo(hypId, gen)
1027 {
1028   _name = "HexaFromSkin_3D";
1029 }
1030
1031 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
1032 {
1033 }
1034
1035 //================================================================================
1036 /*!
1037  * \brief Main method, which generates hexaheda
1038  */
1039 //================================================================================
1040
1041 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1042 {
1043   _Skin skin;
1044   int nbBlocks = skin.findBlocks(aMesh);
1045   if ( nbBlocks == 0 )
1046     return error( skin.error());
1047
1048   vector< vector< const SMDS_MeshNode* > > columns;
1049   int x, xSize, y, ySize, z, zSize;
1050   _Indexer colIndex;
1051
1052   for ( int i = 0; i < nbBlocks; ++i )
1053   {
1054     const _Block& block = skin.getBlock( i );
1055
1056     // ------------------------------------------
1057     // Fill columns of nodes with existing nodes
1058     // ------------------------------------------
1059
1060     xSize = block.getSide(B_BOTTOM).getHoriSize();
1061     ySize = block.getSide(B_BOTTOM).getVertSize();
1062     zSize = block.getSide(B_FRONT ).getVertSize();
1063     int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1064     colIndex = _Indexer( xSize, ySize );
1065     columns.resize( colIndex.size() );
1066
1067     // fill node columns by front and back box sides
1068     for ( x = 0; x < xSize; ++x ) {
1069       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1070       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1071       column0.resize( zSize );
1072       column1.resize( zSize );
1073       for ( z = 0; z < zSize; ++z ) {
1074         column0[ z ] = block.getSide(B_FRONT).node( x, z );
1075         column1[ z ] = block.getSide(B_BACK) .node( x, z );
1076       }
1077     }
1078     // fill node columns by left and right box sides
1079     for ( y = 1; y < ySize-1; ++y ) {
1080       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1081       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1082       column0.resize( zSize );
1083       column1.resize( zSize );
1084       for ( z = 0; z < zSize; ++z ) {
1085         column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1086         column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1087       }
1088     }
1089     // get nodes from top and bottom box sides
1090     for ( x = 1; x < xSize-1; ++x ) {
1091       for ( y = 1; y < ySize-1; ++y ) {
1092         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1093         column.resize( zSize );
1094         column.front() = block.getSide(B_BOTTOM).node( x, y );
1095         column.back()  = block.getSide(B_TOP)   .node( x, y );
1096       }
1097     }
1098
1099     // ----------------------------
1100     // Add internal nodes of a box
1101     // ----------------------------
1102     // projection points of internal nodes on box sub-shapes by which
1103     // coordinates of internal nodes are computed
1104     vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1105
1106     // projections on vertices are constant
1107     pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1108     pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1109     pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1110     pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1111     pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1112     pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1113     pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1114     pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1115
1116     for ( x = 1; x < xSize-1; ++x )
1117     {
1118       gp_XYZ params; // normalized parameters of internal node within a unit box
1119       params.SetCoord( 1, x / double(X) );
1120       for ( y = 1; y < ySize-1; ++y )
1121       {
1122         params.SetCoord( 2, y / double(Y) );
1123         // column to fill during z loop
1124         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1125         // projections on horizontal edges
1126         pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1127         pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1128         pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1129         pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1130         pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1131         pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1132         pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1133         pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1134         // projections on horizontal sides
1135         pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1136         pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP)   .xyz( x, y );
1137         for ( z = 1; z < zSize-1; ++z ) // z loop
1138         {
1139           params.SetCoord( 3, z / double(Z) );
1140           // projections on vertical edges
1141           pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1142           pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1143           pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1144           pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1145           // projections on vertical sides
1146           pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1147           pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1148           pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1149           pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1150
1151           // compute internal node coordinates
1152           gp_XYZ coords;
1153           SMESH_Block::ShellPoint( params, pointOnShape, coords );
1154           column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1155
1156 #ifdef DEB_GRID
1157           // debug
1158           //cout << "----------------------------------------------------------------------"<<endl;
1159           //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1160           //{
1161           //  gp_XYZ p = pointOnShape[ id ];
1162           //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1163           //}
1164           //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1165           //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1166 #endif
1167         }
1168       }
1169     }
1170     // ----------------
1171     // Add hexahedrons
1172     // ----------------
1173
1174     // find out orientation by a least distorted hexahedron (issue 0020855);
1175     // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
1176     double badness = numeric_limits<double>::max();
1177     bool isForw = true;
1178     for ( int xMax = 0; xMax < 2; ++xMax )
1179       for ( int yMax = 0; yMax < 2; ++yMax )
1180         for ( int zMax = 0; zMax < 2; ++zMax )
1181         {
1182           x = xMax ? xSize-1 : 1;
1183           y = yMax ? ySize-1 : 1;
1184           z = zMax ? zSize-1 : 1;
1185           vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
1186           vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x  , y-1 )];
1187           vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y   )];
1188           vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x  , y )];
1189
1190           const SMDS_MeshNode* n000 = col00[z-1];
1191           const SMDS_MeshNode* n100 = col10[z-1];
1192           const SMDS_MeshNode* n010 = col01[z-1];
1193           const SMDS_MeshNode* n110 = col11[z-1];
1194           const SMDS_MeshNode* n001 = col00[z];
1195           const SMDS_MeshNode* n101 = col10[z];
1196           const SMDS_MeshNode* n011 = col01[z];
1197           const SMDS_MeshNode* n111 = col11[z];
1198           SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1199                                           n001,n011,n111,n101);
1200           SMDS_VolumeTool volTool( &probeVolume );
1201           double Nx=0.,Ny=0.,Nz=0.;
1202           for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
1203           {
1204             double nx,ny,nz;
1205             volTool.GetFaceNormal( iFace, nx,ny,nz );
1206             Nx += nx;
1207             Ny += ny;
1208             Nz += nz;
1209           }
1210           double quality = Nx*Nx + Ny*Ny + Nz*Nz;
1211           if ( quality < badness )
1212           {
1213             badness = quality;
1214             isForw = volTool.IsForward();
1215           }
1216         }
1217
1218     // add elements
1219     for ( x = 0; x < xSize-1; ++x ) {
1220       for ( y = 0; y < ySize-1; ++y ) {
1221         vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1222         vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1223         vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1224         vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1225         // bottom face normal of a hexa mush point outside the volume
1226         if ( isForw )
1227           for ( z = 0; z < zSize-1; ++z )
1228             aHelper->AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
1229                                col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1230         else
1231           for ( z = 0; z < zSize-1; ++z )
1232             aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
1233                                col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1234       }
1235     }
1236   } // loop on blocks
1237
1238   return true;
1239 }
1240
1241 //================================================================================
1242 /*!
1243  * \brief Evaluate nb of hexa
1244  */
1245 //================================================================================
1246
1247 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh &         aMesh,
1248                                           const TopoDS_Shape & aShape,
1249                                           MapShapeNbElems&     aResMap)
1250 {
1251   _Skin skin;
1252   int nbBlocks = skin.findBlocks(aMesh);
1253   if ( nbBlocks == 0 )
1254     return error( skin.error());
1255
1256   bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1257
1258   int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1259   vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1260   if ( entity >= (int) nbByType.size() )
1261     nbByType.resize( SMDSEntity_Last, 0 );
1262
1263   for ( int i = 0; i < nbBlocks; ++i )
1264   {
1265     const _Block& block = skin.getBlock( i );
1266
1267     int nbX = block.getSide(B_BOTTOM).getHoriSize();
1268     int nbY = block.getSide(B_BOTTOM).getVertSize();
1269     int nbZ = block.getSide(B_FRONT ).getVertSize();
1270
1271     int nbHexa  = (nbX-1) * (nbY-1) * (nbZ-1);
1272     int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1273     if ( secondOrder )
1274       nbNodes +=
1275         (nbX-2) * (nbY-2) * (nbZ-1) +
1276         (nbX-2) * (nbY-1) * (nbZ-2) +
1277         (nbX-1) * (nbY-2) * (nbZ-2);
1278
1279
1280     nbByType[ entity ] += nbHexa;
1281     nbByType[ SMDSEntity_Node ] += nbNodes;
1282   }
1283
1284   return true;
1285 }
1286
1287 //================================================================================
1288 /*!
1289  * \brief Abstract method must be defined but does nothing
1290  */
1291 //================================================================================
1292
1293 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1294                                                  Hypothesis_Status& aStatus)
1295 {
1296   aStatus = SMESH_Hypothesis::HYP_OK;
1297   return true;
1298 }
1299
1300 //================================================================================
1301 /*!
1302  * \brief Abstract method must be defined but just reports an error as this
1303  *  algo is not intended to work with shapes
1304  */
1305 //================================================================================
1306
1307 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1308 {
1309   return error("Algorithm can't work with geometrical shapes");
1310 }