smesh/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
//
//  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License.
//
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
//
//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
// File      : StdMeshers_HexaFromSkin_3D.cxx
// Created   : Wed Jan 27 12:28:07 2010
// Author    : Edward AGAPOV (eap)


#include "StdMeshers_HexaFromSkin_3D.hxx"

#include "SMDS_VolumeOfNodes.hxx"
#include "SMDS_VolumeTool.hxx"
#include "SMESH_Block.hxx"
#include "SMESH_MesherHelper.hxx"

#include "utilities.h"

// Define error message
#ifdef _DEBUG_
#define BAD_MESH_ERR \
  error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
                      __FILE__ ":" )<<__LINE__)
#else
#define BAD_MESH_ERR \
  error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
#endif

// Debug output
#define _DUMP_(msg) // cout << msg << endl



namespace
{
  enum EBoxSides //!< sides of the block
    {
      B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED
    };
  const char* SBoxSides[] = //!< names of block sides
    {
      "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
    };
  enum EQuadEdge //!< edges of quadrangle side
    {
      Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED
    };

  //================================================================================
  /*!
   * \brief return true if a node is at block corner
   */
  //================================================================================

  bool isCornerNode( const SMDS_MeshNode* n )
  {
    return n && n->NbInverseElements( SMDSAbs_Face ) % 2;
  }

  //================================================================================
  /*!
   * \brief check element type
   */
  //================================================================================

  bool isQuadrangle(const SMDS_MeshElement* e)
  {
    return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 ));
  }

  //================================================================================
  /*!
   * \brief return opposite node of a quadrangle face
   */
  //================================================================================

  const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
  {
    return quad->GetNode( (iNode+2) % 4 );
  }

  //================================================================================
  /*!
   * \brief Convertor of a pair of integers to a sole index
   */
  struct _Indexer
  {
    int _xSize, _ySize;
    _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
    int size() const { return _xSize * _ySize; }
    int operator()(int x, int y) const { return y * _xSize + x; }
  };
  //================================================================================
  /*!
   * \brief Oriented convertor of a pair of integers to a sole index 
   */
  class _OrientedIndexer : public _Indexer
  {
  public:
    enum OriFlags //!< types of block side orientation
      {
        REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
      };
    _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
      _Indexer( indexer._xSize, indexer._ySize ),
      _xSize (indexer._xSize), _ySize(indexer._ySize),
      _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
      _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
      _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
    {
      (*_swapFun)( _xSize, _ySize );
    }
    //!< Return index by XY
    int operator()(int x, int y) const
    {
      (*_xRevFun)( x, const_cast<int&>( _xSize ));
      (*_yRevFun)( y, const_cast<int&>( _ySize ));
      (*_swapFun)( x, y );
      return _Indexer::operator()(x,y);
    }
    //!< Return index for a corner
    int corner(bool xMax, bool yMax) const
    {
      int x = xMax, y = yMax, size = 2;
      (*_xRevFun)( x, size );
      (*_yRevFun)( y, size );
      (*_swapFun)( x, y );
      return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
    }
    int xSize() const { return _xSize; }
    int ySize() const { return _ySize; }
  private:
    _Indexer _indexer;
    int _xSize, _ySize;

    typedef void (*TFun)(int& x, int& y);
    TFun _xRevFun, _yRevFun, _swapFun;
    
    static void lazy(int&, int&) {}
    static void reverse(int& x, int& size) { x = size - x - 1; }
    static void swap(int& x, int& y) { std::swap(x,y); }
  };
  //================================================================================
  /*!
   * \brief Structure corresponding to the meshed side of block
   */
  struct _BlockSide
  {
    vector<const SMDS_MeshNode*> _grid;
    _Indexer                     _index;
    int                          _nbBlocksExpected;
    int                          _nbBlocksFound;

#ifdef _DEBUG_
#define _grid_access_(args) _grid.at( args )
#else
#define _grid_access_(args) _grid[ args ]
#endif
    //!< Return node at XY
    const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); }
    //!< Set node at XY
    void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; }
    //!< Return a corner node
    const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
    {
      return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
    }
    const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
    //!< True if all blocks this side belongs to have beem found
    bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
    //!< Return coordinates of node at XY
    gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); }
    //!< Return gravity center of the four corners and the middle node
    gp_XYZ getGC() const
    {
      gp_XYZ xyz =
        getXYZ( 0, 0 ) +
        getXYZ( _index._xSize-1, 0 ) +
        getXYZ( 0, _index._ySize-1 ) +
        getXYZ( _index._xSize-1, _index._ySize-1 ) +
        getXYZ( _index._xSize/2, _index._ySize/2 );
      return xyz / 5;
    }
    //!< Return number of mesh faces
    int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
  };
  //================================================================================
  /*!
   * \brief _BlockSide with changed orientation
   */
  struct _OrientedBlockSide
  {
    _BlockSide*       _side;
    _OrientedIndexer  _index;

    _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
      _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
    //!< return coordinates by XY
    gp_XYZ xyz(int x, int y) const
    {
      return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) );
    }
    //!< safely return a node by XY
    const SMDS_MeshNode* node(int x, int y) const
    {
      int i = _index( x, y );
      return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
    }
    //!< Return a corner node
    const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
    {
      return _side->_grid_access_( _index.corner( isXMax, isYMax ));
    }
    //!< return its size in nodes
    int getHoriSize() const { return _index.xSize(); }
    int getVertSize() const  { return _index.ySize(); }
    //!< True if _side has been initialized
    operator bool() const { return _side; }
    //! Direct access to _side
    const _BlockSide* operator->() const { return _side; }
    _BlockSide* operator->() { return _side; }
  };
  //================================================================================
  /*!
   * \brief Meshed skin of block
   */
  struct _Block
  {
    _OrientedBlockSide _side[6]; // 6 sides of a sub-block

    const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
  };
  //================================================================================
  /*!
   * \brief Skin mesh possibly containing several meshed blocks
   */
  class _Skin
  {
  public:

    int findBlocks(SMESH_Mesh& mesh);
    //!< return i-th block
    const _Block& getBlock(int i) const { return _blocks[i]; }
    //!< return error description
    const SMESH_Comment& error() const { return _error; }

  private:
    bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad);
    bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
                             const SMDS_MeshNode*    n1,
                             const SMDS_MeshNode*    n2,
                             vector<const SMDS_MeshNode*>& verRow1,
                             vector<const SMDS_MeshNode*>& verRow2,
                             bool alongN1N2 );
    _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
                                      EQuadEdge sharedSideEdge1,
                                      EQuadEdge sharedSideEdge2);
    //!< update own data and data of the side bound to block
    void setSideBoundToBlock( _BlockSide& side )
    {
      side._nbBlocksFound++;
      if ( side.isBound() )
      {
        _corner2sides[ side.getCornerNode(0,0) ].erase( &side );
        _corner2sides[ side.getCornerNode(1,0) ].erase( &side );
        _corner2sides[ side.getCornerNode(0,1) ].erase( &side );
        _corner2sides[ side.getCornerNode(1,1) ].erase( &side );
      }
    }
    //!< store reason of error
    int error(const SMESH_Comment& reason) { _error = reason; return 0; }

    SMESH_Comment      _error;

    list< _BlockSide > _allSides;
    vector< _Block >   _blocks;

    map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
  };

  //================================================================================
  /*!
   * \brief Find and return number of submeshes corresponding to blocks
   */
  //================================================================================

  int _Skin::findBlocks(SMESH_Mesh& mesh)
  {
    SMESHDS_Mesh* meshDS = mesh.GetMeshDS();

    // Find a node at any block corner

    SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
    if ( !nIt->more() ) return error("Empty mesh");

    const SMDS_MeshNode* nCorner = 0;
    while ( nIt->more() )
    {
      nCorner = nIt->next();
      if ( isCornerNode( nCorner ))
        break;
      else
        nCorner = 0;
    }
    if ( !nCorner )
      return BAD_MESH_ERR;

    // --------------------------------------------------------------------
    // Find all block sides starting from mesh faces sharing the corner node
    // --------------------------------------------------------------------

    int nbFacesOnSides = 0;
    TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
    list< const SMDS_MeshNode* > corners( 1, nCorner );
    list< const SMDS_MeshNode* >::iterator corner = corners.begin();
    while ( corner != corners.end() )
    {
      SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
      while ( faceIt->more() )
      {
        const SMDS_MeshElement* face = faceIt->next();
        if ( !cornerFaces.insert( face ).second )
          continue; // already loaded block side

        if ( !isQuadrangle( face ))
          return error("Non-quadrangle elements in the input mesh");

        if ( _allSides.empty() || !_allSides.back()._grid.empty() )
          _allSides.push_back( _BlockSide() );

        _BlockSide& side = _allSides.back();
        if ( !fillSide( side, face ) )
        {
          if ( !_error.empty() )
            return false;
        }
        else
        {
          for ( int isXMax = 0; isXMax < 2; ++isXMax )
            for ( int isYMax = 0; isYMax < 2; ++isYMax )
            {
              const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
              if ( !isCornerNode( nCorner ))
                return BAD_MESH_ERR;
              _corner2sides[ nCorner ].insert( &side );
              corners.push_back( nCorner );
              cornerFaces.insert( side.getCornerFace( nCorner ));
            }
          nbFacesOnSides += side.getNbFaces();
        }
      }
      ++corner;

      // find block sides of other domains if any
      if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
      {
        while ( nIt->more() )
        {
          nCorner = nIt->next();
          if ( isCornerNode( nCorner ))
            corner = corners.insert( corner, nCorner );
        }
        nbFacesOnSides = mesh.NbQuadrangles();
      }
    }
    
    if ( _allSides.empty() )
      return BAD_MESH_ERR;
    if ( _allSides.back()._grid.empty() )
      _allSides.pop_back();

    // ---------------------------
    // Organize sides into blocks
    // ---------------------------

    // analyse sharing of sides by blocks
    int nbBlockSides = 0; // nb of block sides taking into account their sharing
    list < _BlockSide >::iterator sideIt = _allSides.begin();
    for ( ; sideIt != _allSides.end(); ++sideIt )
    {
      _BlockSide& side = *sideIt;
      bool isSharedSide = true;
      for ( int isXMax = 0; isXMax < 2; ++isXMax )
        for ( int isYMax = 0; isYMax < 2; ++isYMax )
          if ( _corner2sides[ side.getCornerNode(isXMax,isYMax) ].size() < 5 )
            isSharedSide = false;
      side._nbBlocksFound = 0;
      side._nbBlocksExpected = isSharedSide ? 2 : 1;
      nbBlockSides += side._nbBlocksExpected;
    }

    // find sides of each block
    int nbBlocks = 0;
    while ( nbBlockSides >= 6 )
    {
      // get any side not bound to all blocks it belongs to
      sideIt = _allSides.begin();
      while ( sideIt != _allSides.end() && sideIt->isBound())
        ++sideIt;

      // start searching for block sides from the got side
      bool ok = true;
      if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
        _blocks.resize( _blocks.size() + 1 );

      _Block& block = _blocks.back();
      block._side[B_FRONT] = &(*sideIt);
      setSideBoundToBlock( *sideIt );
      nbBlockSides--;
      
      // edges of neighbour sides of B_FRONT corresponding to front's edges
      EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
      EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
      for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i )
        ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i]));
      if ( ok )
        ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP ));

      // count the found sides
      _DUMP_(endl);
      for (int i = 0; i < B_UNDEFINED; ++i )
      {
        _DUMP_("** Block side "<< SBoxSides[i]);
        if ( block._side[ i ] )
        {
          if ( ok && i != B_FRONT)
          {
            setSideBoundToBlock( *block._side[ i ]._side );
            nbBlockSides--;
          }
          _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0));
          _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0));    
          _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1));
          _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1));
        }
        else
        {
          _DUMP_("Not found"<<endl);
        }
      }
      if ( !ok )
        block._side[0] = block._side[1] = block._side[2] =
          block._side[3] = block._side[4] = block._side[5] = 0;
      else
        nbBlocks++;
    }
    if ( nbBlocks == 0 && _error.empty() )
      return BAD_MESH_ERR;

    return nbBlocks;
  }

  //================================================================================
  /*!
   * \brief Fill block side data starting from its corner quadrangle
   */
  //================================================================================

  bool _Skin::fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad)
  {
    // Find out size of block side mesured in nodes and by the way find two rows
    // of nodes in two directions.

    int x, y, nbX, nbY;
    const SMDS_MeshElement* firstQuad = cornerQuad;
    {
      // get corner node of cornerQuad
      const SMDS_MeshNode* nCorner = 0;
      for ( int i = 0; i < 4 && !nCorner; ++i )
        if ( isCornerNode( firstQuad->GetNode(i)))
          nCorner = firstQuad->GetNode(i);
      if ( !nCorner ) return false;

      // get a node on block edge
      int iCorner = firstQuad->GetNodeIndex( nCorner );
      const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);

      // find out size of block side
      vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
      if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
           !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
        return false;
      nbX = horRow1.size(), nbY = verRow1.size();

      // store found nodes
      side._index._xSize = horRow1.size();
      side._index._ySize = verRow1.size();
      side._grid.resize( side._index.size(), NULL );

      for ( x = 0; x < horRow1.size(); ++x )
      {
        side.setNode( x, 0, horRow1[x] );
        side.setNode( x, 1, horRow2[x] );
      }
      for ( y = 0; y < verRow1.size(); ++y )
      {
        side.setNode( 0, y, verRow1[y] );
        side.setNode( 1, y, verRow2[y] );
      }
    }
    // Find the rest nodes

    y = 1; // y of the row to fill
    TIDSortedElemSet emptySet, avoidSet;
    while ( ++y < nbY )
    {
      // get next firstQuad in the next row of quadrangles
      //
      //          n2up
      //     o---o               <- y row
      //     |   |
      //     o---o  o  o  o  o   <- found nodes
      //n1down    n2down       
      //
      int i1down, i2down, i2up;
      const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
      const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
      avoidSet.clear(); avoidSet.insert( firstQuad );
      firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
                                                   &i1down, &i2down);
      if ( !isQuadrangle( firstQuad ))
        return BAD_MESH_ERR;

      const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
      avoidSet.clear(); avoidSet.insert( firstQuad );

      // find the rest nodes in the y-th row by faces in the row

      x = 1; 
      while ( ++x < nbX )
      {
        const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
                                                                        avoidSet, &i2up, &i2down);
        if ( !isQuadrangle( quad ))
          return BAD_MESH_ERR;

        n2up   = oppositeNode( quad, i2down );
        n2down = oppositeNode( quad, i2up );
        avoidSet.clear(); avoidSet.insert( quad );

        side.setNode( x, y, n2up );
      }
    }

    return true;
  }

  //================================================================================
  /*!
   * \brief Try to find a block side adjacent to the given side by given edge
   */
  //================================================================================

  _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
                                           EQuadEdge sharedSideEdge1,
                                           EQuadEdge sharedSideEdge2)
  {
    _Block& block = _blocks.back();
    _OrientedBlockSide& side1 = block._side[ startBlockSide ];

    // get corner nodes of the given block edge
    bool xMax1=0, yMax1=0, xMax2=1, yMax2=1;
    switch( sharedSideEdge1 )
    {
    case Q_BOTTOM: yMax2 = 0; break;
    case Q_RIGHT:  xMax1 = 1; break;
    case Q_TOP:    yMax1 = 1; break;
    case Q_LEFT:   xMax2 = 0; break;
    default:
      error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__);
      return 0;
    }
    const SMDS_MeshNode* n1 = side1.cornerNode( xMax1, yMax1);
    const SMDS_MeshNode* n2 = side1.cornerNode( xMax2, yMax2);

    set< _BlockSide* >& sidesWithN1 = _corner2sides[ n1 ];
    set< _BlockSide* >& sidesWithN2 = _corner2sides[ n2 ];
    if ( sidesWithN1.empty() || sidesWithN2.empty() )
    {
      _DUMP_("no sides by nodes "<< n1->GetID() << " " << n2->GetID() );
      return 0;
    }

    // find all sides sharing both nodes n1 and n2
    set< _BlockSide* > sidesOnEdge;
    set< _BlockSide* >::iterator sideIt;
    std::set_intersection( sidesWithN1.begin(), sidesWithN1.end(),
                           sidesWithN2.begin(), sidesWithN2.end(),
                           inserter( sidesOnEdge, sidesOnEdge.begin()));

    // exclude loaded sides of block from sidesOnEdge
    int nbLoadedSides = 0;
    for (int i = 0; i < B_UNDEFINED; ++i )
    {
      if ( block._side[ i ] )
      {
        nbLoadedSides++;
        sidesOnEdge.erase( block._side[ i ]._side );
      }
    }
    int nbSidesOnEdge = sidesOnEdge.size();
    _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge);
    if ( nbSidesOnEdge == 0 )
      return 0;

    _BlockSide* foundSide = 0;
    if ( nbSidesOnEdge == 1 || nbSidesOnEdge == 2 && nbLoadedSides == 1 )
    {
      foundSide = *sidesOnEdge.begin();
    }
    else
    {
      // Select one of found sides most close to startBlockSide

      // gravity center of already loaded block sides
      gp_XYZ gc;
      for (int i = 0; i < B_UNDEFINED; ++i )
        if ( block._side[ i ] )
          gc += block._side[ i ]._side->getGC();
      gc /= nbLoadedSides;

      gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
      gp_Vec p2p1( p1 - p2 );

      const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
      gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
      gp_Vec side1Dir( p1, p1Op );
      
      map < double , _BlockSide* > angleOfSide;
      for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
      {
        _BlockSide* sideI = *sideIt;
        const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
        gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
        gp_Vec sideIDir( p1, p1Op );
        double angle = side1Dir.AngleWithRef( sideIDir, p2p1 );
        if ( angle < 0 ) angle += 2 * PI;
        angleOfSide.insert( make_pair( angle, sideI ));
      }
      gp_Vec gcDir( p1, gc );
      double gcAngle = side1Dir.AngleWithRef( gcDir, p2p1 );
      foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
    }

    // Orient the found side correctly

    // corners of found side corresponding to nodes n1 and n2
    xMax1= yMax1 = 0; xMax2 = yMax2 = 1;
    switch( sharedSideEdge2 )
    {
    case Q_BOTTOM: yMax2 = 0; break;
    case Q_RIGHT:  xMax1 = 1; break;
    case Q_TOP:    yMax1 = 1; break;
    case Q_LEFT:   xMax2 = 0; break;
    default:
      error(SMESH_Comment("Internal error at")<<__FILE__<<":"<<__LINE__);
      return 0;
    }
    for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
    {
      _OrientedBlockSide orientedSide( foundSide, ori );
      const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
      const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
      if ( n1 == n12 && n2 == n22 )
        return orientedSide;
    }
    error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
          << " of side " << startBlockSide
          << " of block " << _blocks.size());
    return 0;
  }

  //================================================================================
  /*!
   * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
   * from the given quadrangle until another block corner encounters.
   *  n1 and n2 are at bottom of quad, n1 is at block corner.
   */
  //================================================================================

  bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement*       quad,
                                  const SMDS_MeshNode*          n1,
                                  const SMDS_MeshNode*          n2,
                                  vector<const SMDS_MeshNode*>& row1,
                                  vector<const SMDS_MeshNode*>& row2,
                                  const bool                    alongN1N2 )
  {
    const SMDS_MeshNode* corner1 = n1;

    // Store nodes of quad in the rows and find new n1 and n2 to get
    // the next face so that new n2 is on block edge
    int i1 = quad->GetNodeIndex( n1 );
    int i2 = quad->GetNodeIndex( n2 );
    row1.clear(); row2.clear();
    row1.push_back( n1 );
    if ( alongN1N2 )
    {
      row1.push_back( n2 );
      row2.push_back( oppositeNode( quad, i2 ));
      row2.push_back( n1 = oppositeNode( quad, i1 ));
    }
    else
    {
      row2.push_back( n2 );
      row1.push_back( n2 = oppositeNode( quad, i2 ));
      row2.push_back( n1 = oppositeNode( quad, i1 ));
    }

    // Find the rest nodes
    TIDSortedElemSet emptySet, avoidSet;
    while ( !isCornerNode( n2 ))
    {
      avoidSet.clear(); avoidSet.insert( quad );
      quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
      if ( !isQuadrangle( quad ))
        return BAD_MESH_ERR;

      row1.push_back( n2 = oppositeNode( quad, i1 ));
      row2.push_back( n1 = oppositeNode( quad, i2 ));
    }
    return n1 != corner1;
  }

  //================================================================================
  /*!
   * \brief Return a corner face by a corner node
   */
  //================================================================================

  const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
  {
    int x, y, isXMax, isYMax, found = 0;
    for ( isXMax = 0; isXMax < 2; ++isXMax )
    {
      for ( isYMax = 0; isYMax < 2; ++isYMax )
      {
        x = isXMax ? _index._xSize-1 : 0;
        y = isYMax ? _index._ySize-1 : 0;
        found = ( getNode(x,y) == cornerNode );
        if ( found ) break;
      }
      if ( found ) break;
    }
    if ( !found ) return 0;
    int dx = isXMax ? -1 : +1;
    int dy = isYMax ? -1 : +1;
    const SMDS_MeshNode* n1 = getNode(x,y);
    const SMDS_MeshNode* n2 = getNode(x+dx,y);
    const SMDS_MeshNode* n3 = getNode(x,y+dy);
    const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
    return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
  }

}

//=======================================================================
//function : StdMeshers_HexaFromSkin_3D
//purpose  : 
//=======================================================================

StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen)
  :SMESH_3D_Algo(hypId, studyId, gen)
{
  MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D");
  _name = "HexaFromSkin_3D";
}

StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
{
  MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D");
}

//================================================================================
/*!
 * \brief Main method, which generates hexaheda
 */
//================================================================================

bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
{
  _Skin skin;
  int nbBlocks = skin.findBlocks(aMesh);
  if ( nbBlocks == 0 )
    return error( skin.error());

  vector< vector< const SMDS_MeshNode* > > columns;
  int x, xSize, y, ySize, z, zSize;
  _Indexer colIndex;

  for ( int i = 0; i < nbBlocks; ++i )
  {
    const _Block& block = skin.getBlock( i );

    // ------------------------------------------
    // Fill columns of nodes with existing nodes
    // ------------------------------------------

    xSize = block.getSide(B_BOTTOM).getHoriSize();
    ySize = block.getSide(B_BOTTOM).getVertSize();
    zSize = block.getSide(B_FRONT ).getVertSize();
    int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
    colIndex = _Indexer( xSize, ySize );
    columns.resize( colIndex.size() );

    // fill node columns by front and back box sides
    for ( x = 0; x < xSize; ++x ) {
      vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
      vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
      column0.resize( zSize );
      column1.resize( zSize );
      for ( z = 0; z < zSize; ++z ) {
        column0[ z ] = block.getSide(B_FRONT).node( x, z );
        column1[ z ] = block.getSide(B_BACK) .node( x, z );
      }
    }
    // fill node columns by left and right box sides
    for ( y = 1; y < ySize-1; ++y ) {
      vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
      vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
      column0.resize( zSize );
      column1.resize( zSize );
      for ( z = 0; z < zSize; ++z ) {
        column0[ z ] = block.getSide(B_LEFT) .node( y, z );
        column1[ z ] = block.getSide(B_RIGHT).node( y, z );
      }
    }
    // get nodes from top and bottom box sides
    for ( x = 1; x < xSize-1; ++x ) {
      for ( y = 1; y < ySize-1; ++y ) {
        vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
        column.resize( zSize );
        column.front() = block.getSide(B_BOTTOM).node( x, y );
        column.back()  = block.getSide(B_TOP)   .node( x, y );
      }
    }

    // ----------------------------
    // Add internal nodes of a box
    // ----------------------------
    // projection points of internal nodes on box subshapes by which
    // coordinates of internal nodes are computed
    vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );

    // projections on vertices are constant
    pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
    pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
    pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
    pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
    pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
    pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
    pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
    pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );

    for ( x = 1; x < xSize-1; ++x )
    {
      gp_XYZ params; // normalized parameters of internal node within a unit box
      params.SetCoord( 1, x / double(X) );
      for ( y = 1; y < ySize-1; ++y )
      {
        params.SetCoord( 2, y / double(Y) );
        // column to fill during z loop
        vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
        // projections on horizontal edges
        pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
        pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
        pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
        pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
        pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
        pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
        pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
        pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
        // projections on horizontal sides
        pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
        pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP)   .xyz( x, y );
        for ( z = 1; z < zSize-1; ++z ) // z loop
        {
          params.SetCoord( 3, z / double(Z) );
          // projections on vertical edges
          pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );    
          pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );    
          pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );    
          pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
          // projections on vertical sides
          pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );    
          pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );    
          pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );    
          pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );

          // compute internal node coordinates
          gp_XYZ coords;
          SMESH_Block::ShellPoint( params, pointOnShape, coords );
          column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );

#ifdef DEB_GRID
          // debug
          //cout << "----------------------------------------------------------------------"<<endl;
          //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
          //{
          //  gp_XYZ p = pointOnShape[ id ];
          //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
          //}
          //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
          //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
#endif
        }
      }
    }
    // ----------------
    // Add hexahedrons
    // ----------------

    // find out orientation
    const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
    const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
    const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
    const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
    const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
    const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
    const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
    const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
    SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
                                    n001,n011,n111,n101);
    bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();

    // add elements
    for ( x = 0; x < xSize-1; ++x ) {
      for ( y = 0; y < ySize-1; ++y ) {
        vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
        vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
        vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
        vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
        // bottom face normal of a hexa mush point outside the volume
        if ( isForw )
          for ( z = 0; z < zSize-1; ++z )
            aHelper->AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
                               col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
        else
          for ( z = 0; z < zSize-1; ++z )
            aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
                               col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
      }
    }
  } // loop on blocks

  return true;
}

//================================================================================
/*!
 * \brief Evaluate nb of hexa
 */
//================================================================================

bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh &         aMesh,
                                          const TopoDS_Shape & aShape,
                                          MapShapeNbElems&     aResMap)
{
  _Skin skin;
  int nbBlocks = skin.findBlocks(aMesh);
  if ( nbBlocks == 0 )
    return error( skin.error());

  bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );

  int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
  vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
  if ( entity >= nbByType.size() )
    nbByType.resize( SMDSEntity_Last, 0 );

  for ( int i = 0; i < nbBlocks; ++i )
  {
    const _Block& block = skin.getBlock( i );

    int nbX = block.getSide(B_BOTTOM).getHoriSize();
    int nbY = block.getSide(B_BOTTOM).getVertSize();
    int nbZ = block.getSide(B_FRONT ).getVertSize();

    int nbHexa  = (nbX-1) * (nbY-1) * (nbZ-1);
    int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
    if ( secondOrder )
      nbNodes +=
        (nbX-2) * (nbY-2) * (nbZ-1) +
        (nbX-2) * (nbY-1) * (nbZ-2) +
        (nbX-1) * (nbY-2) * (nbZ-2);


    nbByType[ entity ] += nbHexa;
    nbByType[ SMDSEntity_Node ] += nbNodes;
  }

  return true;
}

//================================================================================
/*!
 * \brief Abstract method must be defined but does nothing
 */
//================================================================================

bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
                                                 Hypothesis_Status& aStatus)
{
  aStatus = SMESH_Hypothesis::HYP_OK;
  return true;
}

//================================================================================
/*!
 * \brief Abstract method must be defined but just reports an error as this
 *  algo is not intended to work with shapes
 */
//================================================================================

bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
{
  return error("Algorithm can't work with geometrical shapes");
}