src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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