src/StdMeshers/StdMeshers_Hexa_3D.cxx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
  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
  23 //  SMESH SMESH : implementation of SMESH idl descriptions
  24 //  File   : StdMeshers_Hexa_3D.cxx
  25 //           Moved here from SMESH_Hexa_3D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //
  29 #include "StdMeshers_Hexa_3D.hxx"
  30
  31 #include "SMDS_MeshNode.hxx"
  32 #include "SMESH_Comment.hxx"
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_Mesh.hxx"
  35 #include "SMESH_MesherHelper.hxx"
  36 #include "SMESH_subMesh.hxx"
  37 #include "StdMeshers_BlockRenumber.hxx"
  38 #include "StdMeshers_CompositeHexa_3D.hxx"
  39 #include "StdMeshers_FaceSide.hxx"
  40 #include "StdMeshers_HexaFromSkin_3D.hxx"
  41 #include "StdMeshers_Penta_3D.hxx"
  42 #include "StdMeshers_Prism_3D.hxx"
  43 #include "StdMeshers_Quadrangle_2D.hxx"
  44 #include "StdMeshers_ViscousLayers.hxx"
  45
  46 #include <BRep_Tool.hxx>
  47 #include <Bnd_B3d.hxx>
  48 #include <TopExp.hxx>
  49 #include <TopExp_Explorer.hxx>
  50 #include <TopTools_MapOfShape.hxx>
  51 #include <TopTools_SequenceOfShape.hxx>
  52 #include <TopoDS.hxx>
  53
  54 #include "utilities.h"
  55 #include "Utils_ExceptHandlers.hxx"
  56
  57 #include <cstddef>
  58 #include <numeric>
  59
  60 typedef SMESH_Comment TComm;
  61
  62 using namespace std;
  63
  64 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
  65                                                     const TopoDS_Shape &,
  66                                                     SMESH_ProxyMesh* proxyMesh=0);
  67
  68 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
  69                                     MapShapeNbElems &);
  70
  71 //=============================================================================
  72 /*!
  73  * Constructor
  74  */
  75 //=============================================================================
  76
  77 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, SMESH_Gen * gen)
  78   :SMESH_3D_Algo(hypId, gen)
  79 {
  80   _name = "Hexa_3D";
  81   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
  82   _requireShape = false;
  83   _compatibleHypothesis.push_back("ViscousLayers");
  84   _compatibleHypothesis.push_back("BlockRenumber");
  85   _quadAlgo = new StdMeshers_Quadrangle_2D( gen->GetANewId(), _gen );
  86 }
  87
  88 //=============================================================================
  89 /*!
  90  * Destructor
  91  */
  92 //=============================================================================
  93
  94 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
  95 {
  96   delete _quadAlgo;
  97   _quadAlgo = 0;
  98 }
  99
 100 //=============================================================================
 101 /*!
 102  * Retrieves defined hypotheses
 103  */
 104 //=============================================================================
 105
 106 bool StdMeshers_Hexa_3D::CheckHypothesis
 107                          (SMESH_Mesh&                          aMesh,
 108                           const TopoDS_Shape&                  aShape,
 109                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 110 {
 111   // check nb of faces in the shape
 112 /*  PAL16229
 113   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
 114   int nbFaces = 0;
 115   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
 116     if ( ++nbFaces > 6 )
 117       break;
 118   if ( nbFaces != 6 )
 119     return false;
 120 */
 121
 122   _viscousLayersHyp = nullptr;
 123   _blockRenumberHyp = nullptr;
 124
 125   const list<const SMESHDS_Hypothesis*>& hyps =
 126     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 127   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 128   if ( h == hyps.end())
 129   {
 130     aStatus = SMESH_Hypothesis::HYP_OK;
 131     return true;
 132   }
 133
 134   // only StdMeshers_ViscousLayers can be used
 135   aStatus = HYP_OK;
 136   for ( ; h != hyps.end(); ++h )
 137   {
 138     if ( !_viscousLayersHyp &&
 139          (_viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h )))
 140       continue;
 141     if ( !_blockRenumberHyp &&
 142          (_blockRenumberHyp = dynamic_cast< const StdMeshers_BlockRenumber*> ( *h )))
 143       continue;
 144     break;
 145   }
 146   if ((int) hyps.size() != (bool)_viscousLayersHyp + (bool)_blockRenumberHyp )
 147     aStatus = HYP_INCOMPATIBLE;
 148   else
 149   {
 150     if ( _viscousLayersHyp )
 151       if ( !error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus )))
 152         aStatus = HYP_BAD_PARAMETER;
 153
 154     if ( _blockRenumberHyp && aStatus == HYP_OK )
 155       error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape ));
 156   }
 157
 158   return aStatus == HYP_OK;
 159 }
 160
 161 namespace
 162 {
 163   //=============================================================================
 164
 165   typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
 166   typedef std::vector<gp_XYZ>                 TXYZColumn;
 167
 168   // symbolic names of box sides
 169   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
 170
 171   // symbolic names of sides of quadrangle
 172   enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
 173
 174   enum EAxes{ COO_X=1, COO_Y, COO_Z };
 175
 176   //=============================================================================
 177   /*!
 178    * \brief Container of nodes of structured mesh on a qudrangular geom FACE
 179    */
 180   struct _FaceGrid
 181   {
 182     // face sides
 183     FaceQuadStructPtr _quad;
 184
 185     // map of (node parameter on EDGE) to (column (vector) of nodes)
 186     TParam2ColumnMap _u2nodesMap;
 187
 188     // node column's taken from _u2nodesMap taking into account sub-shape orientation
 189     vector<TNodeColumn> _columns;
 190
 191     // columns of normalized parameters of nodes within the unitary cube
 192     vector<TXYZColumn> _ijkColumns;
 193
 194     // geometry of a cube side
 195     TopoDS_Face _sideF;
 196
 197     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
 198     {
 199       return _columns[iCol][iRow];
 200     }
 201     gp_XYZ GetXYZ(int iCol, int iRow) const
 202     {
 203       return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
 204     }
 205     gp_XYZ& GetIJK(int iCol, int iRow)
 206     {
 207       return _ijkColumns[iCol][iRow];
 208     }
 209   };
 210
 211   //================================================================================
 212   /*!
 213    * \brief Converter of a pair of integers to a sole index
 214    */
 215   struct _Indexer
 216   {
 217     int _xSize, _ySize;
 218     _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
 219     int size() const { return _xSize * _ySize; }
 220     int operator()(const int x, const int y) const { return y * _xSize + x; }
 221   };
 222
 223   //================================================================================
 224   /*!
 225    * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
 226    *  the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
 227    */
 228   FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSidePtr side,
 229                                        FaceQuadStructPtr      quad[ 6 ])
 230   {
 231     FaceQuadStructPtr foundQuad;
 232     for ( int i = 1; i < 6; ++i )
 233     {
 234       if ( !quad[i] ) continue;
 235       for ( size_t iS = 0; iS < quad[i]->side.size(); ++iS )
 236       {
 237         const StdMeshers_FaceSidePtr side2 = quad[i]->side[iS];
 238         if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
 239               side->FirstVertex().IsSame( side2->LastVertex() ))
 240             &&
 241             ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
 242               side->LastVertex().IsSame( side2->LastVertex() ))
 243             )
 244         {
 245           if ( iS != Q_BOTTOM )
 246           {
 247             vector< FaceQuadStruct::Side > newSides;
 248             for ( size_t j = iS; j < quad[i]->side.size(); ++j )
 249               newSides.push_back( quad[i]->side[j] );
 250             for ( size_t j = 0; j < iS; ++j )
 251               newSides.push_back( quad[i]->side[j] );
 252             quad[i]->side.swap( newSides );
 253           }
 254           foundQuad.swap(quad[i]);
 255           return foundQuad;
 256         }
 257       }
 258     }
 259     return foundQuad;
 260   }
 261
 262   //================================================================================
 263   /*!
 264    * \brief Put quads to aCubeSide in the order of enum EBoxSides
 265    */
 266   //================================================================================
 267
 268   bool arrangeQuads( FaceQuadStructPtr quad[ 6 ], _FaceGrid aCubeSide[ 6 ], bool reverseBottom )
 269   {
 270     swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
 271     if ( reverseBottom )
 272       swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the bottom normal inside cube
 273             aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
 274
 275     aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
 276     aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
 277     aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP   ], quad );
 278     aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT  ], quad );
 279     if ( aCubeSide[B_FRONT ]._quad )
 280       aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
 281
 282     for ( int i = 1; i < 6; ++i )
 283       if ( !aCubeSide[i]._quad )
 284         return false;
 285     return true;
 286   }
 287
 288   //================================================================================
 289   /*!
 290    * \brief Rearrange block sides according to StdMeshers_BlockRenumber hypothesis
 291    */
 292   //================================================================================
 293
 294   bool arrangeForRenumber( _FaceGrid      blockSide[ 6 ],
 295                            TopoDS_Vertex& v000,
 296                            TopoDS_Vertex& v001 )
 297   {
 298     std::swap( blockSide[B_BOTTOM]._quad->side[ Q_RIGHT],// restore after arrangeQuads()
 299                blockSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
 300
 301     // find v000
 302     TopTools_MapOfShape cornerVertices;
 303     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_BOTTOM].grid->LastVertex()  );
 304     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_BOTTOM].grid->FirstVertex() );
 305     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_TOP   ].grid->LastVertex()  );
 306     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_TOP   ].grid->FirstVertex() );
 307     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_BOTTOM].grid->FirstVertex() );
 308     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_BOTTOM].grid->LastVertex()  );
 309     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_TOP   ].grid->FirstVertex() );
 310     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_TOP   ].grid->LastVertex()  );
 311
 312     if ( v000.IsNull() )
 313     {
 314       // block CS is not defined;
 315       // renumber only if the block has an edge parallel to an axis of global CS
 316
 317       v000 = StdMeshers_RenumberHelper::GetVertex000( cornerVertices );
 318     }
 319
 320     Bnd_B3d bbox;
 321     for ( auto it = cornerVertices.cbegin(); it != cornerVertices.cend(); ++it )
 322       bbox.Add( BRep_Tool::Pnt( TopoDS::Vertex( *it )));
 323     double tol = 1e-5 * Sqrt( bbox.SquareExtent() );
 324
 325     // get block edges starting at v000
 326
 327     std::vector< StdMeshers_FaceSidePtr > edgesAtV000;
 328     std::vector< gp_Vec >                 edgeDir;
 329     std::vector< int >                    iParallel; // 0 - none, 1 - X, 2 - Y, 3 - Z
 330     TopTools_MapOfShape                   lastVertices;
 331     for ( int iQ = 0; iQ < 6; ++iQ )
 332     {
 333       FaceQuadStructPtr quad = blockSide[iQ]._quad;
 334       for ( size_t iS = 0; iS < quad->side.size() &&  edgesAtV000.size() < 3; ++iS )
 335       {
 336         StdMeshers_FaceSidePtr edge = quad->side[iS];
 337         TopoDS_Vertex v1 = edge->FirstVertex(), v2 = edge->LastVertex();
 338         if (( v1.IsSame( v000 ) && !lastVertices.Contains( v2 )) ||
 339             ( v2.IsSame( v000 ) && !lastVertices.Contains( v1 )))
 340         {
 341           bool reverse = v2.IsSame( v000 );
 342           if ( reverse )
 343             std::swap( v1, v2 );
 344           lastVertices.Add( v2 );
 345
 346           edgesAtV000.push_back( edge );
 347
 348           gp_Pnt pf = BRep_Tool::Pnt( v1 );
 349           gp_Pnt pl = BRep_Tool::Pnt( v2 );
 350           gp_Vec vec( pf, pl );
 351           edgeDir.push_back( vec );
 352
 353           iParallel.push_back( 0 );
 354           if ( !v001.IsNull() )
 355           {
 356             if ( v001.IsSame( v2 ))
 357               iParallel.back() = 3;
 358           }
 359           else
 360           {
 361             bool isStraight = true;
 362             for ( int iE = 0; iE < edge->NbEdges() &&  isStraight; ++iE )
 363               isStraight = SMESH_Algo::IsStraight( edge->Edge( iE ));
 364
 365             // is parallel to a GCS axis?
 366             if ( isStraight )
 367             {
 368               int nbDiff = (( Abs( vec.X() ) > tol ) +
 369                             ( Abs( vec.Y() ) > tol ) +
 370                             ( Abs( vec.Z() ) > tol ) );
 371               if ( nbDiff == 1 )
 372                 iParallel.back() = ( Abs( vec.X() ) > tol ) ? 1 : ( Abs( vec.Y() ) > tol ) ? 2 : 3;
 373             }
 374             else
 375             {
 376               edgeDir.back() = gp_Vec( pf, edge->Value3d( reverse ? 0.99 : 0.01 ));
 377             }
 378           }
 379         }
 380       }
 381     }
 382     if ( std::accumulate( iParallel.begin(), iParallel.end(), 0 ) == 0 )
 383       return false;
 384
 385     // find edge OZ and edge OX
 386     StdMeshers_FaceSidePtr edgeOZ, edgeOX;
 387     auto iZIt = std::find( iParallel.begin(), iParallel.end(), 3 );
 388     if ( iZIt != iParallel.end() )
 389     {
 390       int i = std::distance( iParallel.begin(), iZIt );
 391       edgeOZ = edgesAtV000[ i ];
 392       int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() );
 393       int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() );
 394       if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 )
 395         std::swap( iE1, iE2 );
 396       edgeOX = edgesAtV000[ iE1 ];
 397     }
 398     else
 399     {
 400       for ( size_t i = 0; i < edgesAtV000.size(); ++i )
 401       {
 402         if ( !iParallel[ i ] )
 403           continue;
 404         int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() );
 405         int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() );
 406         if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 )
 407           std::swap( iE1, iE2 );
 408         edgeOZ = edgesAtV000[ iParallel[i] == 1 ? iE2 : iE1 ];
 409         edgeOX = edgesAtV000[ iParallel[i] == 1 ? i : iE1 ];
 410         break;
 411       }
 412     }
 413
 414     if ( !edgeOZ || !edgeOX )
 415       return false;
 416
 417     TopoDS_Vertex v100 = edgeOX->LastVertex();
 418     if ( v100.IsSame( v000 ))
 419       v100 = edgeOX->FirstVertex();
 420
 421     // Find the left quad, one including v000 but not v100
 422
 423     for ( int iQ = 0; iQ < 6; ++iQ )
 424     {
 425       FaceQuadStructPtr quad = blockSide[iQ]._quad;
 426       bool hasV000 = false, hasV100 = false;
 427       for ( size_t iS = 0; iS < quad->side.size(); ++iS )
 428       {
 429         StdMeshers_FaceSidePtr edge = quad->side[iS];
 430         if ( edge->FirstVertex().IsSame( v000 ) || edge->LastVertex().IsSame( v000 ))
 431           hasV000 = true;
 432         if ( edge->FirstVertex().IsSame( v100 ) || edge->LastVertex().IsSame( v100 ))
 433           hasV100 = true;
 434       }
 435       if ( hasV000 && !hasV100 )
 436       {
 437         // orient the left quad
 438         for ( int i = 0; i < 4; ++i )
 439         {
 440           if ( quad->side[Q_BOTTOM].grid->Edge(0).IsSame( edgeOZ->Edge(0) ))
 441             break;
 442           quad->shift( 1, true );
 443         }
 444
 445         FaceQuadStructPtr quads[ 6 ];
 446         quads[0].swap( blockSide[iQ]._quad );
 447         for ( int i = 1, j = 0; i < 6; ++i, ++j )
 448           if ( blockSide[ j ]._quad )
 449             quads[ i ].swap( blockSide[ j ]._quad );
 450           else
 451             --i;
 452
 453         return arrangeQuads( quads, blockSide, false/* true*/ );
 454       }
 455     }
 456     return false;
 457   }
 458
 459   //================================================================================
 460   /*!
 461    * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
 462    */
 463   //================================================================================
 464
 465   bool beginsAtSide( const _FaceGrid&     sideGrid1,
 466                      const _FaceGrid&     sideGrid2,
 467                      SMESH_ProxyMesh::Ptr proxymesh )
 468   {
 469     const TNodeColumn& col0  = sideGrid2._u2nodesMap.begin()->second;
 470     const TNodeColumn& col1  = sideGrid2._u2nodesMap.rbegin()->second;
 471     const SMDS_MeshNode* n00 = col0.front();
 472     const SMDS_MeshNode* n01 = col0.back();
 473     const SMDS_MeshNode* n10 = col1.front();
 474     const SMDS_MeshNode* n11 = col1.back();
 475     const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
 476     if ( proxymesh )
 477     {
 478       n00 = proxymesh->GetProxyNode( n00 );
 479       n10 = proxymesh->GetProxyNode( n10 );
 480       n01 = proxymesh->GetProxyNode( n01 );
 481       n11 = proxymesh->GetProxyNode( n11 );
 482       n   = proxymesh->GetProxyNode( n );
 483     }
 484     return ( n == n00 || n == n01 || n == n10 || n == n11 );
 485   }
 486
 487   //================================================================================
 488   /*!
 489    * \brief Fill in _FaceGrid::_ijkColumns
 490    *  \param [in,out] fg - a _FaceGrid
 491    *  \param [in] i1 - coordinate index along _columns
 492    *  \param [in] i2 - coordinate index along _columns[i]
 493    *  \param [in] v3 - value of the constant parameter
 494    */
 495   //================================================================================
 496
 497   void computeIJK( _FaceGrid& fg, int i1, int i2, double v3 )
 498   {
 499     gp_XYZ ijk( v3, v3, v3 );
 500     const size_t nbCol = fg._columns.size();
 501     const size_t nbRow = fg._columns[0].size();
 502
 503     fg._ijkColumns.resize( nbCol );
 504     for ( size_t i = 0; i < nbCol; ++i )
 505       fg._ijkColumns[ i ].resize( nbRow, ijk );
 506
 507     vector< double > len( nbRow );
 508     len[0] = 0;
 509     for ( size_t i = 0; i < nbCol; ++i )
 510     {
 511       gp_Pnt pPrev = fg.GetXYZ( i, 0 );
 512       for ( size_t j = 1; j < nbRow; ++j )
 513       {
 514         gp_Pnt p = fg.GetXYZ( i, j );
 515         len[ j ] = len[ j-1 ] + p.Distance( pPrev );
 516         pPrev = p;
 517       }
 518       for ( size_t j = 0; j < nbRow; ++j )
 519         fg.GetIJK( i, j ).SetCoord( i2, len[ j ]/len.back() );
 520     }
 521
 522     len.resize( nbCol );
 523     for ( size_t j = 0; j < nbRow; ++j )
 524     {
 525       gp_Pnt pPrev = fg.GetXYZ( 0, j );
 526       for ( size_t i = 1; i < nbCol; ++i )
 527       {
 528         gp_Pnt p = fg.GetXYZ( i, j );
 529         len[ i ] = len[ i-1 ] + p.Distance( pPrev );
 530         pPrev = p;
 531       }
 532       for ( size_t i = 0; i < nbCol; ++i )
 533         fg.GetIJK( i, j ).SetCoord( i1, len[ i ]/len.back() );
 534     }
 535   }
 536 }
 537
 538 //=============================================================================
 539 /*!
 540  * Generates hexahedron mesh on hexaedron like form using algorithm from
 541  * "Application de l'interpolation transfinie � la cr�ation de maillages
 542  *  C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
 543  *  et hexaedres d�form�s."
 544  * Alain PERONNET - 8 janvier 1999
 545  */
 546 //=============================================================================
 547
 548 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
 549                                  const TopoDS_Shape & aShape)
 550 {
 551   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 552   //Unexpect aCatch(SalomeException);
 553   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 554
 555   // Shape verification
 556   // ----------------------
 557
 558   // shape must be a solid (or a shell) with 6 faces
 559   TopExp_Explorer exp(aShape,TopAbs_SHELL);
 560   if ( !exp.More() )
 561     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
 562   if ( exp.Next(), exp.More() )
 563     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
 564
 565   TopTools_IndexedMapOfShape FF, EE;
 566   TopExp::MapShapes( aShape, TopAbs_FACE, FF);
 567   if ( FF.Extent() != 6)
 568   {
 569     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), _gen);
 570     compositeHexa.SetHypothesis( _blockRenumberHyp );
 571     if ( !compositeHexa.Compute( aMesh, aShape ))
 572       return error( compositeHexa.GetComputeError() );
 573
 574     return _blockRenumberHyp ? error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape )) : true;
 575   }
 576
 577   // Find sides of a cube
 578   // ---------------------
 579
 580   // tool creating quadratic elements if needed
 581   SMESH_MesherHelper helper (aMesh);
 582   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 583
 584   TopExp::MapShapes( aShape, TopAbs_EDGE, EE );
 585   SMESH_MesherHelper* faceHelper = ( EE.Size() == 12 ) ? 0 : &helper;
 586
 587   FaceQuadStructPtr quad[ 6 ];
 588   for ( int i = 0; i < 6; ++i )
 589   {
 590     if ( faceHelper )
 591       faceHelper->SetSubShape( FF( i+1 ));
 592     if ( !( quad[i] = FaceQuadStructPtr( _quadAlgo->CheckNbEdges( aMesh, FF( i+1 ),
 593                                                                   /*considerMesh=*/true,
 594                                                                   faceHelper))))
 595       return error( _quadAlgo->GetComputeError() );
 596     if ( quad[i]->side.size() != 4 )
 597       return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
 598   }
 599
 600   // put quads in a proper order
 601   _FaceGrid aCubeSide[ 6 ];
 602   if ( !arrangeQuads( quad, aCubeSide, true ))
 603     return error( COMPERR_BAD_SHAPE );
 604
 605
 606   // Make viscous layers
 607   // --------------------
 608
 609   SMESH_ProxyMesh::Ptr proxymesh;
 610   if ( _viscousLayersHyp )
 611   {
 612     proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
 613     if ( !proxymesh )
 614       return false;
 615   }
 616
 617   // Check if there are triangles on cube sides
 618   // -------------------------------------------
 619
 620   if ( aMesh.NbTriangles() > 0 )
 621   {
 622     for ( int i = 0; i < 6; ++i )
 623     {
 624       const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
 625       const SMESHDS_SubMesh* smDS =
 626         proxymesh ? proxymesh->GetSubMesh( sideF ) : meshDS->MeshElements( sideF );
 627       if ( !SMESH_MesherHelper::IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE,
 628                                                     /*nullSubMeshRes=*/false ))
 629       {
 630         SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 631         return error( err );
 632       }
 633     }
 634   }
 635
 636   // Arrange sides according to _blockRenumberHyp
 637   bool toRenumber = _blockRenumberHyp;
 638   if ( toRenumber )
 639   {
 640     TopoDS_Vertex v000, v001;
 641     _blockRenumberHyp->IsSolidIncluded( aMesh, aShape, v000, v001 );
 642
 643     toRenumber = arrangeForRenumber( aCubeSide, v000, v001 );
 644
 645     if ( toRenumber )
 646     {
 647       meshDS->Modified();
 648       meshDS->CompactMesh(); // remove numbering holes
 649     }
 650   }
 651
 652   // Check presence of regular grid mesh on FACEs of the cube
 653   // ------------------------------------------------------------
 654
 655   for ( int i = 0; i < 6; ++i )
 656   {
 657     const TopoDS_Face& F = aCubeSide[i]._quad->face;
 658     StdMeshers_FaceSidePtr baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
 659     list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
 660
 661     // assure correctness of node positions on baseE:
 662     // helper.GetNodeU() will fix positions if they are wrong
 663     helper.ToFixNodeParameters( true );
 664     for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
 665     {
 666       const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
 667       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
 668       {
 669         bool ok;
 670         helper.SetSubShape( baseE );
 671         SMDS_ElemIteratorPtr eIt = smDS->GetElements();
 672         while ( eIt->more() )
 673         {
 674           const SMDS_MeshElement* e = eIt->next();
 675           // expect problems on a composite side
 676           try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
 677           catch (...) {}
 678           try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
 679           catch (...) {}
 680         }
 681       }
 682     }
 683
 684     // load grid
 685     bool ok =
 686       helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
 687     if ( ok )
 688     {
 689       // check if the loaded grid corresponds to nb of quadrangles on the FACE
 690       const SMESHDS_SubMesh* faceSubMesh =
 691         proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
 692       const smIdType nbQuads = faceSubMesh->NbElements();
 693       const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
 694       const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
 695       ok = ( nbQuads == nbHor * nbVer );
 696     }
 697     if ( !ok )
 698     {
 699       SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 700       return error( err );
 701     }
 702   }
 703
 704   // Orient loaded grids of cube sides along axis of the unitary cube coord system
 705   bool isReverse[6];
 706   isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
 707   isReverse[B_TOP   ] = beginsAtSide( aCubeSide[B_TOP   ], aCubeSide[B_RIGHT ], proxymesh );
 708   isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
 709   isReverse[B_BACK  ] = beginsAtSide( aCubeSide[B_BACK  ], aCubeSide[B_RIGHT ], proxymesh );
 710   isReverse[B_LEFT  ] = beginsAtSide( aCubeSide[B_LEFT  ], aCubeSide[B_BACK  ], proxymesh );
 711   isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK  ], proxymesh );
 712   for ( int i = 0; i < 6; ++i )
 713   {
 714     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
 715
 716     size_t iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
 717     size_t*  pi = isReverse[i] ? &iRev : &iFwd;
 718     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
 719     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
 720       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
 721
 722     aCubeSide[i]._u2nodesMap.clear();
 723   }
 724
 725   if ( proxymesh )
 726     for ( int i = 0; i < 6; ++i )
 727       for ( size_t j = 0; j < aCubeSide[i]._columns.size(); ++j)
 728         for ( size_t k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
 729         {
 730           const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
 731           n = proxymesh->GetProxyNode( n );
 732         }
 733
 734   // 4) Create internal nodes of the cube
 735   // -------------------------------------
 736
 737   helper.SetSubShape( aShape );
 738   helper.SetElementsOnShape(true);
 739
 740   // shortcuts to sides
 741   _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
 742   _FaceGrid* fRight  = & aCubeSide[ B_RIGHT  ];
 743   _FaceGrid* fTop    = & aCubeSide[ B_TOP    ];
 744   _FaceGrid* fLeft   = & aCubeSide[ B_LEFT   ];
 745   _FaceGrid* fFront  = & aCubeSide[ B_FRONT  ];
 746   _FaceGrid* fBack   = & aCubeSide[ B_BACK   ];
 747
 748   // cube size measured in nb of nodes
 749   size_t x, xSize = fBottom->_columns.size() , X = xSize - 1;
 750   size_t y, ySize = fLeft->_columns.size()   , Y = ySize - 1;
 751   size_t z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
 752
 753   // check sharing of FACEs (IPAL54417)
 754   if ( fFront ->_columns.size()    != xSize ||
 755        fBack  ->_columns.size()    != xSize ||
 756        fTop   ->_columns.size()    != xSize ||
 757
 758        fRight ->_columns.size()    != ySize ||
 759        fTop   ->_columns[0].size() != ySize ||
 760        fBottom->_columns[0].size() != ySize ||
 761
 762        fRight ->_columns[0].size() != zSize ||
 763        fFront ->_columns[0].size() != zSize ||
 764        fBack  ->_columns[0].size() != zSize )
 765     return error( COMPERR_BAD_SHAPE, "Not sewed faces" );
 766
 767   // columns of internal nodes "rising" from nodes of fBottom
 768   _Indexer colIndex( xSize, ySize );
 769   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
 770
 771   // fill node columns by front and back box sides
 772   for ( x = 0; x < xSize; ++x ) {
 773     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
 774     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
 775     column0.resize( zSize );
 776     column1.resize( zSize );
 777     for ( z = 0; z < zSize; ++z ) {
 778       column0[ z ] = fFront->GetNode( x, z );
 779       column1[ z ] = fBack ->GetNode( x, z );
 780     }
 781   }
 782   // fill node columns by left and right box sides
 783   for ( y = 1; y < ySize-1; ++y ) {
 784     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
 785     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
 786     column0.resize( zSize );
 787     column1.resize( zSize );
 788     for ( z = 0; z < zSize; ++z ) {
 789       column0[ z ] = fLeft ->GetNode( y, z );
 790       column1[ z ] = fRight->GetNode( y, z );
 791     }
 792   }
 793   // get nodes from top and bottom box sides
 794   for ( x = 1; x < xSize-1; ++x ) {
 795     for ( y = 1; y < ySize-1; ++y ) {
 796       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 797       column.resize( zSize );
 798       column.front() = fBottom->GetNode( x, y );
 799       column.back()  = fTop   ->GetNode( x, y );
 800     }
 801   }
 802
 803   // compute normalized parameters of nodes on sides (PAL23189)
 804   computeIJK( *fBottom, COO_X, COO_Y, /*z=*/0. );
 805   computeIJK( *fRight,  COO_Y, COO_Z, /*x=*/1. );
 806   computeIJK( *fTop,    COO_X, COO_Y, /*z=*/1. );
 807   computeIJK( *fLeft,   COO_Y, COO_Z, /*x=*/0. );
 808   computeIJK( *fFront,  COO_X, COO_Z, /*y=*/0. );
 809   computeIJK( *fBack,   COO_X, COO_Z, /*y=*/1. );
 810
 811   StdMeshers_RenumberHelper renumHelper( aMesh, _blockRenumberHyp );
 812
 813   // projection points of the internal node on cube sub-shapes by which
 814   // coordinates of the internal node are computed
 815   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
 816
 817   // projections on vertices are constant
 818   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
 819   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
 820   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
 821   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
 822   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
 823   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
 824   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
 825   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
 826
 827   gp_XYZ params; // normalized parameters of an internal node within the unit box
 828
 829   if ( toRenumber )
 830     for ( y = 0; y < ySize; ++y )
 831     {
 832       vector< const SMDS_MeshNode* >& column0y = columns[ colIndex( 0, y )];
 833       for ( z = 0; z < zSize; ++z )
 834         renumHelper.AddReplacingNode( column0y[ z ] );
 835     }
 836
 837   for ( x = 1; x < xSize-1; ++x )
 838   {
 839     if ( toRenumber )
 840     {
 841       vector< const SMDS_MeshNode* >& columnX0 = columns[ colIndex( x, 0 )];
 842       for ( z = 0; z < zSize; ++z )
 843         renumHelper.AddReplacingNode( columnX0[ z ] );
 844     }
 845
 846     const double rX = x / double(X);
 847     for ( y = 1; y < ySize-1; ++y )
 848     {
 849       const double rY = y / double(Y);
 850
 851       // a column to fill in during z loop
 852       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 853       // projection points on horizontal edges
 854       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
 855       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
 856       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
 857       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
 858       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
 859       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
 860       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
 861       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
 862       // projection points on horizontal faces
 863       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
 864       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
 865
 866       if ( toRenumber )
 867         renumHelper.AddReplacingNode( column[ 0 ] );
 868
 869       for ( z = 1; z < zSize-1; ++z ) // z loop
 870       {
 871         const double rZ = z / double(Z);
 872
 873         const gp_XYZ& pBo = fBottom->GetIJK( x, y );
 874         const gp_XYZ& pTo = fTop   ->GetIJK( x, y );
 875         const gp_XYZ& pFr = fFront ->GetIJK( x, z );
 876         const gp_XYZ& pBa = fBack  ->GetIJK( x, z );
 877         const gp_XYZ& pLe = fLeft  ->GetIJK( y, z );
 878         const gp_XYZ& pRi = fRight ->GetIJK( y, z );
 879         params.SetCoord( 1, 0.5 * ( pBo.X() * ( 1. - rZ ) + pTo.X() * rZ  +
 880                                     pFr.X() * ( 1. - rY ) + pBa.X() * rY ));
 881         params.SetCoord( 2, 0.5 * ( pBo.Y() * ( 1. - rZ ) + pTo.Y() * rZ  +
 882                                     pLe.Y() * ( 1. - rX ) + pRi.Y() * rX ));
 883         params.SetCoord( 3, 0.5 * ( pFr.Z() * ( 1. - rY ) + pBa.Z() * rY  +
 884                                     pLe.Z() * ( 1. - rX ) + pRi.Z() * rX ));
 885
 886         // projection points on vertical edges
 887         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
 888         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
 889         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
 890         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
 891         // projection points on vertical faces
 892         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
 893         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
 894         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
 895         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
 896
 897         // compute internal node coordinates
 898         gp_XYZ coords;
 899         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
 900         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
 901
 902       } // z loop
 903       if ( toRenumber )
 904         renumHelper.AddReplacingNode( column[ Z ] );
 905
 906     } // y loop
 907     if ( toRenumber )
 908     {
 909       vector< const SMDS_MeshNode* >& columnX0 = columns[ colIndex( x, Y )];
 910       for ( z = 0; z < zSize; ++z )
 911         renumHelper.AddReplacingNode( columnX0[ z ] );
 912     }
 913   } // x loop
 914
 915   if ( toRenumber )
 916     for ( y = 0; y < ySize; ++y )
 917     {
 918       vector< const SMDS_MeshNode* >& columnXy = columns[ colIndex( X, y )];
 919       for ( z = 0; z < zSize; ++z )
 920         renumHelper.AddReplacingNode( columnXy[ z ] );
 921     }
 922
 923   // side data no more needed, free memory
 924   for ( int i = 0; i < 6; ++i )
 925     SMESHUtils::FreeVector( aCubeSide[i]._columns );
 926
 927   // 5) Create hexahedrons
 928   // ---------------------
 929
 930   for ( x = 0; x < xSize-1; ++x ) {
 931     for ( y = 0; y < ySize-1; ++y ) {
 932       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
 933       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
 934       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
 935       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
 936       for ( z = 0; z < zSize-1; ++z )
 937       {
 938         // bottom face normal of a hexa mush point outside the volume
 939         if ( toRenumber )
 940           helper.AddVolume(col00[z], col01[z], col01[z+1], col00[z+1],
 941                            col10[z], col11[z], col11[z+1], col10[z+1]);
 942         else
 943           helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
 944                            col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
 945       }
 946     }
 947   }
 948
 949   if ( toRenumber )
 950     renumHelper.DoReplaceNodes();
 951
 952
 953   if ( _blockRenumberHyp )
 954   {
 955     return error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape ));
 956   }
 957
 958   return true;
 959 }
 960
 961 //=============================================================================
 962 /*!
 963  *  Evaluate
 964  */
 965 //=============================================================================
 966
 967 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
 968                                   const TopoDS_Shape & aShape,
 969                                   MapShapeNbElems& aResMap)
 970 {
 971   vector < SMESH_subMesh * >meshFaces;
 972   TopTools_SequenceOfShape aFaces;
 973   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 974     aFaces.Append(exp.Current());
 975     SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 976     ASSERT(aSubMesh);
 977     meshFaces.push_back(aSubMesh);
 978   }
 979   if (meshFaces.size() != 6) {
 980     //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
 981     static StdMeshers_CompositeHexa_3D compositeHexa(-10, aMesh.GetGen());
 982     return compositeHexa.Evaluate(aMesh, aShape, aResMap);
 983   }
 984
 985   int i = 0;
 986   for(; i<6; i++) {
 987     //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
 988     TopoDS_Shape aFace = aFaces.Value(i+1);
 989     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
 990     if( !algo ) {
 991       std::vector<smIdType> aResVec(SMDSEntity_Last);
 992       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 993       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 994       aResMap.insert(std::make_pair(sm,aResVec));
 995       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 996       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 997       return false;
 998     }
 999     string algoName = algo->GetName();
1000     bool isAllQuad = false;
1001     if (algoName == "Quadrangle_2D") {
1002       MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
1003       if( anIt == aResMap.end() ) continue;
1004       std::vector<smIdType> aVec = (*anIt).second;
1005       smIdType nbtri = std::max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1006       if( nbtri == 0 )
1007         isAllQuad = true;
1008     }
1009     if ( ! isAllQuad ) {
1010       return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
1011     }
1012   }
1013
1014   // find number of 1d elems for 1 face
1015   int nb1d = 0;
1016   TopTools_MapOfShape Edges1;
1017   bool IsQuadratic = false;
1018   bool IsFirst = true;
1019   for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
1020     Edges1.Add(exp.Current());
1021     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
1022     if( sm ) {
1023       MapShapeNbElemsItr anIt = aResMap.find(sm);
1024       if( anIt == aResMap.end() ) continue;
1025       std::vector<smIdType> aVec = (*anIt).second;
1026       nb1d += std::max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1027       if(IsFirst) {
1028         IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1029         IsFirst = false;
1030       }
1031     }
1032   }
1033   // find face opposite to 1 face
1034   int OppNum = 0;
1035   for(i=2; i<=6; i++) {
1036     bool IsOpposite = true;
1037     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
1038       if( Edges1.Contains(exp.Current()) ) {
1039         IsOpposite = false;
1040         break;
1041       }
1042     }
1043     if(IsOpposite) {
1044       OppNum = i;
1045       break;
1046     }
1047   }
1048   // find number of 2d elems on side faces
1049   int nb2d = 0;
1050   for(i=2; i<=6; i++) {
1051     if( i == OppNum ) continue;
1052     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
1053     if( anIt == aResMap.end() ) continue;
1054     std::vector<smIdType> aVec = (*anIt).second;
1055     nb2d += std::max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1056   }
1057
1058   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
1059   std::vector<smIdType> aVec = (*anIt).second;
1060   smIdType nb2d_face0 = std::max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1061   smIdType nb0d_face0 = aVec[SMDSEntity_Node];
1062
1063   std::vector<smIdType> aResVec(SMDSEntity_Last);
1064   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1065   if(IsQuadratic) {
1066     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
1067     smIdType nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
1068     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
1069   }
1070   else {
1071     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
1072     aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
1073   }
1074   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1075   aResMap.insert(std::make_pair(sm,aResVec));
1076
1077   return true;
1078 }
1079
1080 //================================================================================
1081 /*!
1082  * \brief Computes hexahedral mesh from 2D mesh of block
1083  */
1084 //================================================================================
1085
1086 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1087 {
1088   static StdMeshers_HexaFromSkin_3D * algo = 0;
1089   if ( !algo ) {
1090     SMESH_Gen* gen = aMesh.GetGen();
1091     algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), gen );
1092   }
1093   algo->InitComputeError();
1094   algo->Compute( aMesh, aHelper );
1095   return error( algo->GetComputeError());
1096 }
1097
1098 //================================================================================
1099 /*!
1100  * \brief Return true if the algorithm can mesh this shape
1101  *  \param [in] aShape - shape to check
1102  *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
1103  *              else, returns OK if at least one shape is OK
1104  */
1105 //================================================================================
1106
1107 bool StdMeshers_Hexa_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1108 {
1109   TopExp_Explorer exp0( aShape, TopAbs_SOLID );
1110   if ( !exp0.More() ) return false;
1111
1112   for ( ; exp0.More(); exp0.Next() )
1113   {
1114     int nbFoundShells = 0;
1115     TopExp_Explorer exp1( exp0.Current(), TopAbs_SHELL );
1116     for ( ; exp1.More(); exp1.Next(), ++nbFoundShells)
1117       if ( nbFoundShells == 2 ) break;
1118     if ( nbFoundShells != 1 ) {
1119       if ( toCheckAll ) return false;
1120       continue;
1121     }
1122     exp1.Init( exp0.Current(), TopAbs_FACE );
1123     int nbEdges = SMESH_MesherHelper::Count( exp1.Current(), TopAbs_EDGE, /*ignoreSame=*/true );
1124     bool ok = ( nbEdges > 3 );
1125     if ( toCheckAll && !ok ) return false;
1126     if ( !toCheckAll && ok ) return true;
1127   }
1128   return toCheckAll;
1129 }
1130
1131 //=======================================================================
1132 //function : ComputePentahedralMesh
1133 //purpose  :
1134 //=======================================================================
1135
1136 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &          aMesh,
1137                                              const TopoDS_Shape &  aShape,
1138                                              SMESH_ProxyMesh*      proxyMesh)
1139 {
1140   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
1141   if ( proxyMesh )
1142   {
1143     err->myName = COMPERR_BAD_INPUT_MESH;
1144     err->myComment = "Can't build pentahedral mesh on viscous layers";
1145     return err;
1146   }
1147   bool bOK;
1148   StdMeshers_Penta_3D anAlgo;
1149   //
1150   bOK=anAlgo.Compute(aMesh, aShape);
1151   //
1152   err = anAlgo.GetComputeError();
1153   //
1154   if ( !bOK && anAlgo.ErrorStatus() == 5 )
1155   {
1156     static StdMeshers_Prism_3D * aPrism3D = 0;
1157     if ( !aPrism3D ) {
1158       SMESH_Gen* gen = aMesh.GetGen();
1159       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
1160     }
1161     SMESH_Hypothesis::Hypothesis_Status aStatus;
1162     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
1163       aPrism3D->InitComputeError();
1164       bOK = aPrism3D->Compute( aMesh, aShape );
1165       err = aPrism3D->GetComputeError();
1166     }
1167   }
1168   return err;
1169 }
1170
1171
1172 //=======================================================================
1173 //function : EvaluatePentahedralMesh
1174 //purpose  :
1175 //=======================================================================
1176
1177 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
1178                              const TopoDS_Shape & aShape,
1179                              MapShapeNbElems& aResMap)
1180 {
1181   StdMeshers_Penta_3D anAlgo;
1182   bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
1183
1184   //err = anAlgo.GetComputeError();
1185   //if ( !bOK && anAlgo.ErrorStatus() == 5 )
1186   if( !bOK ) {
1187     static StdMeshers_Prism_3D * aPrism3D = 0;
1188     if ( !aPrism3D ) {
1189       SMESH_Gen* gen = aMesh.GetGen();
1190       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
1191     }
1192     SMESH_Hypothesis::Hypothesis_Status aStatus;
1193     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
1194       return aPrism3D->Evaluate(aMesh, aShape, aResMap);
1195     }
1196   }
1197
1198   return bOK;
1199 }