src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx

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