src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx

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