src/StdMeshers/StdMeshers_Hexa_3D.cxx

   1 // Copyright (C) 2007-2014  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 : implementaion 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 "StdMeshers_CompositeHexa_3D.hxx"
  32 #include "StdMeshers_FaceSide.hxx"
  33 #include "StdMeshers_HexaFromSkin_3D.hxx"
  34 #include "StdMeshers_Penta_3D.hxx"
  35 #include "StdMeshers_Prism_3D.hxx"
  36 #include "StdMeshers_Quadrangle_2D.hxx"
  37 #include "StdMeshers_ViscousLayers.hxx"
  38
  39 #include "SMESH_Comment.hxx"
  40 #include "SMESH_Gen.hxx"
  41 #include "SMESH_Mesh.hxx"
  42 #include "SMESH_MesherHelper.hxx"
  43 #include "SMESH_subMesh.hxx"
  44
  45 #include "SMDS_MeshNode.hxx"
  46
  47 #include <TopExp.hxx>
  48 #include <TopExp_Explorer.hxx>
  49 #include <TopTools_SequenceOfShape.hxx>
  50 #include <TopTools_MapOfShape.hxx>
  51 #include <TopoDS.hxx>
  52
  53 #include "utilities.h"
  54 #include "Utils_ExceptHandlers.hxx"
  55
  56 typedef SMESH_Comment TComm;
  57
  58 using namespace std;
  59
  60 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
  61                                                     const TopoDS_Shape &,
  62                                                     SMESH_ProxyMesh* proxyMesh=0);
  63
  64 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
  65                                     MapShapeNbElems &);
  66
  67 //=============================================================================
  68 /*!
  69  * Constructor
  70  */
  71 //=============================================================================
  72
  73 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
  74   :SMESH_3D_Algo(hypId, studyId, gen)
  75 {
  76   MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
  77   _name = "Hexa_3D";
  78   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
  79   _requireShape = false;
  80   _compatibleHypothesis.push_back("ViscousLayers");
  81 }
  82
  83 //=============================================================================
  84 /*!
  85  * Destructor
  86  */
  87 //=============================================================================
  88
  89 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
  90 {
  91   MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D");
  92 }
  93
  94 //=============================================================================
  95 /*!
  96  * Retrieves defined hypotheses
  97  */
  98 //=============================================================================
  99
 100 bool StdMeshers_Hexa_3D::CheckHypothesis
 101                          (SMESH_Mesh&                          aMesh,
 102                           const TopoDS_Shape&                  aShape,
 103                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 104 {
 105   // check nb of faces in the shape
 106 /*  PAL16229
 107   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
 108   int nbFaces = 0;
 109   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
 110     if ( ++nbFaces > 6 )
 111       break;
 112   if ( nbFaces != 6 )
 113     return false;
 114 */
 115
 116   _viscousLayersHyp = NULL;
 117
 118   const list<const SMESHDS_Hypothesis*>& hyps =
 119     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 120   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 121   if ( h == hyps.end())
 122   {
 123     aStatus = SMESH_Hypothesis::HYP_OK;
 124     return true;
 125   }
 126
 127   aStatus = HYP_OK;
 128   for ( ; h != hyps.end(); ++h )
 129   {
 130     string hypName = (*h)->GetName();
 131     if ( find( _compatibleHypothesis.begin(),_compatibleHypothesis.end(),hypName )
 132          != _compatibleHypothesis.end() )
 133     {
 134       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 135     }
 136     else
 137     {
 138       aStatus = HYP_INCOMPATIBLE;
 139     }
 140   }
 141
 142   if ( !_viscousLayersHyp )
 143     aStatus = HYP_INCOMPATIBLE;
 144
 145   return aStatus == HYP_OK;
 146 }
 147
 148 namespace
 149 {
 150   //=============================================================================
 151
 152   typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
 153
 154   // symbolic names of box sides
 155   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
 156
 157   // symbolic names of sides of quadrangle
 158   enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
 159
 160   //=============================================================================
 161   /*!
 162    * \brief Container of nodes of structured mesh on a qudrangular geom FACE
 163    */
 164   struct _FaceGrid
 165   {
 166     // face sides
 167     FaceQuadStructPtr _quad;
 168
 169     // map of (node parameter on EDGE) to (column (vector) of nodes)
 170     TParam2ColumnMap _u2nodesMap;
 171
 172     // node column's taken form _u2nodesMap taking into account sub-shape orientation
 173     vector<TNodeColumn> _columns;
 174
 175     // geometry of a cube side
 176     TopoDS_Face _sideF;
 177
 178     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
 179     {
 180       return _columns[iCol][iRow];
 181     }
 182     gp_XYZ GetXYZ(int iCol, int iRow) const
 183     {
 184       return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
 185     }
 186   };
 187
 188   //================================================================================
 189   /*!
 190    * \brief Convertor of a pair of integers to a sole index
 191    */
 192   struct _Indexer
 193   {
 194     int _xSize, _ySize;
 195     _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
 196     int size() const { return _xSize * _ySize; }
 197     int operator()(const int x, const int y) const { return y * _xSize + x; }
 198   };
 199
 200   //================================================================================
 201   /*!
 202    * \brief Appends a range of node columns from a map to another map
 203    */
 204   template< class TMapIterator >
 205   void append( TParam2ColumnMap& toMap, TMapIterator from, TMapIterator to )
 206   {
 207     const SMDS_MeshNode* lastNode = toMap.rbegin()->second[0];
 208     const SMDS_MeshNode* firstNode = from->second[0];
 209     if ( lastNode == firstNode )
 210       from++;
 211     double u = toMap.rbegin()->first;
 212     for (; from != to; ++from )
 213     {
 214       u += 1;
 215       TParam2ColumnMap::iterator u2nn = toMap.insert( toMap.end(), make_pair ( u, TNodeColumn()));
 216       u2nn->second.swap( from->second );
 217     }
 218   }
 219
 220   //================================================================================
 221   /*!
 222    * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
 223    *  the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
 224    */
 225   FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSidePtr side,
 226                                        FaceQuadStructPtr      quad[ 6 ])
 227   {
 228     FaceQuadStructPtr foundQuad;
 229     for ( int i = 1; i < 6; ++i )
 230     {
 231       if ( !quad[i] ) continue;
 232       for ( unsigned iS = 0; iS < quad[i]->side.size(); ++iS )
 233       {
 234         const StdMeshers_FaceSidePtr side2 = quad[i]->side[iS];
 235         if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
 236               side->FirstVertex().IsSame( side2->LastVertex() ))
 237             &&
 238             ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
 239               side->LastVertex().IsSame( side2->LastVertex() ))
 240             )
 241         {
 242           if ( iS != Q_BOTTOM )
 243           {
 244             vector< FaceQuadStruct::Side > newSides;
 245             for ( unsigned j = iS; j < quad[i]->side.size(); ++j )
 246               newSides.push_back( quad[i]->side[j] );
 247             for ( unsigned j = 0; j < iS; ++j )
 248               newSides.push_back( quad[i]->side[j] );
 249             quad[i]->side.swap( newSides );
 250           }
 251           foundQuad.swap(quad[i]);
 252           return foundQuad;
 253         }
 254       }
 255     }
 256     return foundQuad;
 257   }
 258   //================================================================================
 259   /*!
 260    * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
 261    */
 262   //================================================================================
 263
 264   bool beginsAtSide( const _FaceGrid&     sideGrid1,
 265                      const _FaceGrid&     sideGrid2,
 266                      SMESH_ProxyMesh::Ptr proxymesh )
 267   {
 268     const TNodeColumn& col0  = sideGrid2._u2nodesMap.begin()->second;
 269     const TNodeColumn& col1  = sideGrid2._u2nodesMap.rbegin()->second;
 270     const SMDS_MeshNode* n00 = col0.front();
 271     const SMDS_MeshNode* n01 = col0.back();
 272     const SMDS_MeshNode* n10 = col1.front();
 273     const SMDS_MeshNode* n11 = col1.back();
 274     const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
 275     if ( proxymesh )
 276     {
 277       n00 = proxymesh->GetProxyNode( n00 );
 278       n10 = proxymesh->GetProxyNode( n10 );
 279       n01 = proxymesh->GetProxyNode( n01 );
 280       n11 = proxymesh->GetProxyNode( n11 );
 281       n   = proxymesh->GetProxyNode( n );
 282     }
 283     return ( n == n00 || n == n01 || n == n10 || n == n11 );
 284   }
 285 }
 286
 287 //=============================================================================
 288 /*!
 289  * Generates hexahedron mesh on hexaedron like form using algorithm from
 290  * "Application de l'interpolation transfinie � la cr�ation de maillages
 291  *  C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
 292  *  et hexaedres d�form�s."
 293  * Alain PERONNET - 8 janvier 1999
 294  */
 295 //=============================================================================
 296
 297 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
 298                                  const TopoDS_Shape & aShape)
 299 {
 300   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 301   //Unexpect aCatch(SalomeException);
 302   MESSAGE("StdMeshers_Hexa_3D::Compute");
 303   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 304
 305   // Shape verification
 306   // ----------------------
 307
 308   // shape must be a solid (or a shell) with 6 faces
 309   TopExp_Explorer exp(aShape,TopAbs_SHELL);
 310   if ( !exp.More() )
 311     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
 312   if ( exp.Next(), exp.More() )
 313     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
 314
 315   TopTools_IndexedMapOfShape FF;
 316   TopExp::MapShapes( aShape, TopAbs_FACE, FF);
 317   if ( FF.Extent() != 6)
 318   {
 319     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen);
 320     if ( !compositeHexa.Compute( aMesh, aShape ))
 321       return error( compositeHexa.GetComputeError() );
 322     return true;
 323   }
 324
 325   // Find sides of a cube
 326   // ---------------------
 327
 328   FaceQuadStructPtr quad[ 6 ];
 329   StdMeshers_Quadrangle_2D quadAlgo( _gen->GetANewId(), GetStudyId(), _gen);
 330   for ( int i = 0; i < 6; ++i )
 331   {
 332     if ( !( quad[i] = FaceQuadStructPtr( quadAlgo.CheckNbEdges( aMesh, FF( i+1 )))))
 333       return error( quadAlgo.GetComputeError() );
 334     if ( quad[i]->side.size() != 4 )
 335       return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
 336   }
 337
 338   _FaceGrid aCubeSide[ 6 ];
 339
 340   swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
 341   swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the normal of bottom quad inside cube
 342         aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
 343
 344   aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
 345   aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
 346   aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP   ], quad );
 347   aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT  ], quad );
 348   if ( aCubeSide[B_FRONT ]._quad )
 349     aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
 350
 351   for ( int i = 1; i < 6; ++i )
 352     if ( !aCubeSide[i]._quad )
 353       return error( COMPERR_BAD_SHAPE );
 354
 355   // Make viscous layers
 356   // --------------------
 357
 358   SMESH_ProxyMesh::Ptr proxymesh;
 359   if ( _viscousLayersHyp )
 360   {
 361     proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
 362     if ( !proxymesh )
 363       return false;
 364   }
 365
 366   // Check if there are triangles on cube sides
 367   // -------------------------------------------
 368
 369   if ( aMesh.NbTriangles() > 0 )
 370   {
 371     for ( int i = 0; i < 6; ++i )
 372     {
 373       const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
 374       const SMESHDS_SubMesh* smDS =
 375         proxymesh ? proxymesh->GetSubMesh( sideF ) : meshDS->MeshElements( sideF );
 376       if ( !SMESH_MesherHelper::IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE,
 377                                                     /*nullSubMeshRes=*/false ))
 378       {
 379         SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 380         return error( err );
 381       }
 382     }
 383   }
 384
 385   // Check presence of regular grid mesh on FACEs of the cube
 386   // ------------------------------------------------------------
 387
 388   // tool creating quadratic elements if needed
 389   SMESH_MesherHelper helper (aMesh);
 390   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 391
 392   for ( int i = 0; i < 6; ++i )
 393   {
 394     const TopoDS_Face& F = aCubeSide[i]._quad->face;
 395     StdMeshers_FaceSidePtr baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
 396     list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
 397
 398     // assure correctness of node positions on baseE:
 399     // helper.GetNodeU() will fix positions if they are wrong
 400     helper.ToFixNodeParameters( true );
 401     for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
 402     {
 403       const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
 404       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
 405       {
 406         bool ok;
 407         helper.SetSubShape( baseE );
 408         SMDS_ElemIteratorPtr eIt = smDS->GetElements();
 409         while ( eIt->more() )
 410         {
 411           const SMDS_MeshElement* e = eIt->next();
 412           // expect problems on a composite side
 413           try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
 414           catch (...) {}
 415           try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
 416           catch (...) {}
 417         }
 418       }
 419     }
 420
 421     // load grid
 422     bool ok =
 423       helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
 424     if ( ok )
 425     {
 426       // check if the loaded grid corresponds to nb of quadrangles on the FACE
 427       const SMESHDS_SubMesh* faceSubMesh =
 428         proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
 429       const int nbQuads = faceSubMesh->NbElements();
 430       const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
 431       const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
 432       ok = ( nbQuads == nbHor * nbVer );
 433     }
 434     if ( !ok )
 435     {
 436       SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 437       return error( err );
 438     }
 439   }
 440
 441   // Orient loaded grids of cube sides along axis of the unitary cube coord system
 442   bool isReverse[6];
 443   isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
 444   isReverse[B_TOP   ] = beginsAtSide( aCubeSide[B_TOP   ], aCubeSide[B_RIGHT ], proxymesh );
 445   isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
 446   isReverse[B_BACK  ] = beginsAtSide( aCubeSide[B_BACK  ], aCubeSide[B_RIGHT ], proxymesh );
 447   isReverse[B_LEFT  ] = beginsAtSide( aCubeSide[B_LEFT  ], aCubeSide[B_BACK  ], proxymesh );
 448   isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK  ], proxymesh );
 449   for ( int i = 0; i < 6; ++i )
 450   {
 451     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
 452
 453     int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
 454     int* pi = isReverse[i] ? &iRev : &iFwd;
 455     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
 456     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
 457       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
 458
 459     aCubeSide[i]._u2nodesMap.clear();
 460   }
 461
 462   if ( proxymesh )
 463     for ( int i = 0; i < 6; ++i )
 464       for ( unsigned j = 0; j < aCubeSide[i]._columns.size(); ++j)
 465         for ( unsigned k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
 466         {
 467           const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
 468           n = proxymesh->GetProxyNode( n );
 469         }
 470
 471   // 4) Create internal nodes of the cube
 472   // -------------------------------------
 473
 474   helper.SetSubShape( aShape );
 475   helper.SetElementsOnShape(true);
 476
 477   // shortcuts to sides
 478   _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
 479   _FaceGrid* fRight  = & aCubeSide[ B_RIGHT  ];
 480   _FaceGrid* fTop    = & aCubeSide[ B_TOP    ];
 481   _FaceGrid* fLeft   = & aCubeSide[ B_LEFT   ];
 482   _FaceGrid* fFront  = & aCubeSide[ B_FRONT  ];
 483   _FaceGrid* fBack   = & aCubeSide[ B_BACK   ];
 484
 485   // cube size measured in nb of nodes
 486   int x, xSize = fBottom->_columns.size() , X = xSize - 1;
 487   int y, ySize = fLeft->_columns.size()   , Y = ySize - 1;
 488   int z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
 489
 490   // columns of internal nodes "rising" from nodes of fBottom
 491   _Indexer colIndex( xSize, ySize );
 492   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
 493
 494   // fill node columns by front and back box sides
 495   for ( x = 0; x < xSize; ++x ) {
 496     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
 497     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
 498     column0.resize( zSize );
 499     column1.resize( zSize );
 500     for ( z = 0; z < zSize; ++z ) {
 501       column0[ z ] = fFront->GetNode( x, z );
 502       column1[ z ] = fBack ->GetNode( x, z );
 503     }
 504   }
 505   // fill node columns by left and right box sides
 506   for ( y = 1; y < ySize-1; ++y ) {
 507     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
 508     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
 509     column0.resize( zSize );
 510     column1.resize( zSize );
 511     for ( z = 0; z < zSize; ++z ) {
 512       column0[ z ] = fLeft ->GetNode( y, z );
 513       column1[ z ] = fRight->GetNode( y, z );
 514     }
 515   }
 516   // get nodes from top and bottom box sides
 517   for ( x = 1; x < xSize-1; ++x ) {
 518     for ( y = 1; y < ySize-1; ++y ) {
 519       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 520       column.resize( zSize );
 521       column.front() = fBottom->GetNode( x, y );
 522       column.back()  = fTop   ->GetNode( x, y );
 523     }
 524   }
 525
 526   // projection points of the internal node on cube sub-shapes by which
 527   // coordinates of the internal node are computed
 528   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
 529
 530   // projections on vertices are constant
 531   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
 532   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
 533   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
 534   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
 535   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
 536   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
 537   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
 538   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
 539
 540   for ( x = 1; x < xSize-1; ++x )
 541   {
 542     gp_XYZ params; // normalized parameters of internal node within a unit box
 543     params.SetCoord( 1, x / double(X) );
 544     for ( y = 1; y < ySize-1; ++y )
 545     {
 546       params.SetCoord( 2, y / double(Y) );
 547       // a column to fill in during z loop
 548       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 549       // projection points on horizontal edges
 550       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
 551       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
 552       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
 553       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
 554       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
 555       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
 556       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
 557       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
 558       // projection points on horizontal faces
 559       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
 560       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
 561       for ( z = 1; z < zSize-1; ++z ) // z loop
 562       {
 563         params.SetCoord( 3, z / double(Z) );
 564         // projection points on vertical edges
 565         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
 566         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
 567         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
 568         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
 569         // projection points on vertical faces
 570         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
 571         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
 572         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
 573         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
 574
 575         // compute internal node coordinates
 576         gp_XYZ coords;
 577         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
 578         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
 579
 580       }
 581     }
 582   }
 583
 584   // side data no more needed, free memory
 585   for ( int i = 0; i < 6; ++i )
 586     aCubeSide[i]._columns.clear();
 587
 588   // 5) Create hexahedrons
 589   // ---------------------
 590
 591   for ( x = 0; x < xSize-1; ++x ) {
 592     for ( y = 0; y < ySize-1; ++y ) {
 593       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
 594       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
 595       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
 596       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
 597       for ( z = 0; z < zSize-1; ++z )
 598       {
 599         // bottom face normal of a hexa mush point outside the volume
 600         helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
 601                          col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
 602       }
 603     }
 604   }
 605   return true;
 606 }
 607
 608 //=============================================================================
 609 /*!
 610  *  Evaluate
 611  */
 612 //=============================================================================
 613
 614 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
 615                                   const TopoDS_Shape & aShape,
 616                                   MapShapeNbElems& aResMap)
 617 {
 618   vector < SMESH_subMesh * >meshFaces;
 619   TopTools_SequenceOfShape aFaces;
 620   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 621     aFaces.Append(exp.Current());
 622     SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 623     ASSERT(aSubMesh);
 624     meshFaces.push_back(aSubMesh);
 625   }
 626   if (meshFaces.size() != 6) {
 627     //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
 628     static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
 629     return compositeHexa.Evaluate(aMesh, aShape, aResMap);
 630   }
 631
 632   int i = 0;
 633   for(; i<6; i++) {
 634     //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
 635     TopoDS_Shape aFace = aFaces.Value(i+1);
 636     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
 637     if( !algo ) {
 638       std::vector<int> aResVec(SMDSEntity_Last);
 639       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 640       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 641       aResMap.insert(std::make_pair(sm,aResVec));
 642       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 643       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 644       return false;
 645     }
 646     string algoName = algo->GetName();
 647     bool isAllQuad = false;
 648     if (algoName == "Quadrangle_2D") {
 649       MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
 650       if( anIt == aResMap.end() ) continue;
 651       std::vector<int> aVec = (*anIt).second;
 652       int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 653       if( nbtri == 0 )
 654         isAllQuad = true;
 655     }
 656     if ( ! isAllQuad ) {
 657       return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
 658     }
 659   }
 660
 661   // find number of 1d elems for 1 face
 662   int nb1d = 0;
 663   TopTools_MapOfShape Edges1;
 664   bool IsQuadratic = false;
 665   bool IsFirst = true;
 666   for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
 667     Edges1.Add(exp.Current());
 668     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
 669     if( sm ) {
 670       MapShapeNbElemsItr anIt = aResMap.find(sm);
 671       if( anIt == aResMap.end() ) continue;
 672       std::vector<int> aVec = (*anIt).second;
 673       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 674       if(IsFirst) {
 675         IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 676         IsFirst = false;
 677       }
 678     }
 679   }
 680   // find face opposite to 1 face
 681   int OppNum = 0;
 682   for(i=2; i<=6; i++) {
 683     bool IsOpposite = true;
 684     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
 685       if( Edges1.Contains(exp.Current()) ) {
 686         IsOpposite = false;
 687         break;
 688       }
 689     }
 690     if(IsOpposite) {
 691       OppNum = i;
 692       break;
 693     }
 694   }
 695   // find number of 2d elems on side faces
 696   int nb2d = 0;
 697   for(i=2; i<=6; i++) {
 698     if( i == OppNum ) continue;
 699     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
 700     if( anIt == aResMap.end() ) continue;
 701     std::vector<int> aVec = (*anIt).second;
 702     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 703   }
 704
 705   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
 706   std::vector<int> aVec = (*anIt).second;
 707   int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 708   int nb0d_face0 = aVec[SMDSEntity_Node];
 709
 710   std::vector<int> aResVec(SMDSEntity_Last);
 711   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 712   if(IsQuadratic) {
 713     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
 714     int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
 715     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
 716   }
 717   else {
 718     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
 719     aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
 720   }
 721   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 722   aResMap.insert(std::make_pair(sm,aResVec));
 723
 724   return true;
 725 }
 726
 727 //================================================================================
 728 /*!
 729  * \brief Computes hexahedral mesh from 2D mesh of block
 730  */
 731 //================================================================================
 732
 733 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
 734 {
 735   static StdMeshers_HexaFromSkin_3D * algo = 0;
 736   if ( !algo ) {
 737     SMESH_Gen* gen = aMesh.GetGen();
 738     algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen );
 739   }
 740   algo->InitComputeError();
 741   algo->Compute( aMesh, aHelper );
 742   return error( algo->GetComputeError());
 743 }
 744
 745 //================================================================================
 746 /*!
 747  * \brief Return true if the algorithm can mesh this shape
 748  *  \param [in] aShape - shape to check
 749  *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
 750  *              else, returns OK if at least one shape is OK
 751  */
 752 //================================================================================
 753
 754 bool StdMeshers_Hexa_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
 755 {
 756   TopExp_Explorer exp0( aShape, TopAbs_SOLID );
 757   if ( !exp0.More() ) return false;
 758
 759   for ( ; exp0.More(); exp0.Next() )
 760   {
 761     int nbFoundShells = 0;
 762     TopExp_Explorer exp1( exp0.Current(), TopAbs_SHELL );
 763     for ( ; exp1.More(); exp1.Next(), ++nbFoundShells)
 764       if ( nbFoundShells == 2 ) break;
 765     if ( nbFoundShells != 1 ) {
 766       if ( toCheckAll ) return false;
 767       continue;
 768     }
 769     exp1.Init( exp0.Current(), TopAbs_FACE );
 770     int nbEdges = SMESH_MesherHelper::Count( exp1.Current(), TopAbs_EDGE, /*ignoreSame=*/true );
 771     bool ok = ( nbEdges > 3 );
 772     if ( toCheckAll && !ok ) return false;
 773     if ( !toCheckAll && ok ) return true;
 774   }
 775   return toCheckAll;
 776 };
 777
 778 //=======================================================================
 779 //function : ComputePentahedralMesh
 780 //purpose  :
 781 //=======================================================================
 782
 783 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &          aMesh,
 784                                              const TopoDS_Shape &  aShape,
 785                                              SMESH_ProxyMesh*      proxyMesh)
 786 {
 787   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
 788   if ( proxyMesh )
 789   {
 790     err->myName = COMPERR_BAD_INPUT_MESH;
 791     err->myComment = "Can't build pentahedral mesh on viscous layers";
 792     return err;
 793   }
 794   bool bOK;
 795   StdMeshers_Penta_3D anAlgo;
 796   //
 797   bOK=anAlgo.Compute(aMesh, aShape);
 798   //
 799   err = anAlgo.GetComputeError();
 800   //
 801   if ( !bOK && anAlgo.ErrorStatus() == 5 )
 802   {
 803     static StdMeshers_Prism_3D * aPrism3D = 0;
 804     if ( !aPrism3D ) {
 805       SMESH_Gen* gen = aMesh.GetGen();
 806       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
 807     }
 808     SMESH_Hypothesis::Hypothesis_Status aStatus;
 809     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
 810       aPrism3D->InitComputeError();
 811       bOK = aPrism3D->Compute( aMesh, aShape );
 812       err = aPrism3D->GetComputeError();
 813     }
 814   }
 815   return err;
 816 }
 817
 818
 819 //=======================================================================
 820 //function : EvaluatePentahedralMesh
 821 //purpose  :
 822 //=======================================================================
 823
 824 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
 825                              const TopoDS_Shape & aShape,
 826                              MapShapeNbElems& aResMap)
 827 {
 828   StdMeshers_Penta_3D anAlgo;
 829   bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
 830
 831   //err = anAlgo.GetComputeError();
 832   //if ( !bOK && anAlgo.ErrorStatus() == 5 )
 833   if( !bOK ) {
 834     static StdMeshers_Prism_3D * aPrism3D = 0;
 835     if ( !aPrism3D ) {
 836       SMESH_Gen* gen = aMesh.GetGen();
 837       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
 838     }
 839     SMESH_Hypothesis::Hypothesis_Status aStatus;
 840     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
 841       return aPrism3D->Evaluate(aMesh, aShape, aResMap);
 842     }
 843   }
 844
 845   return bOK;
 846 }