src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx

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