src/StdMeshers/StdMeshers_Hexa_3D.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 //  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  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_MeshElement.hxx"
  46 #include "SMDS_MeshNode.hxx"
  47 #include "SMDS_FacePosition.hxx"
  48 #include "SMDS_VolumeTool.hxx"
  49 #include "SMDS_VolumeOfNodes.hxx"
  50
  51 #include <TopExp.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  54 #include <TopTools_ListIteratorOfListOfShape.hxx>
  55 #include <TopTools_ListOfShape.hxx>
  56 #include <TopTools_SequenceOfShape.hxx>
  57 #include <TopTools_MapOfShape.hxx>
  58 #include <TopoDS.hxx>
  59 #include <gp_Pnt2d.hxx>
  60
  61 #include "utilities.h"
  62 #include "Utils_ExceptHandlers.hxx"
  63
  64 typedef SMESH_Comment TComm;
  65
  66 using namespace std;
  67
  68 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
  69                                                     const TopoDS_Shape &,
  70                                                     SMESH_ProxyMesh* proxyMesh=0);
  71
  72 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
  73                                     MapShapeNbElems &);
  74
  75 //=============================================================================
  76 /*!
  77  * Constructor
  78  */
  79 //=============================================================================
  80
  81 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
  82   :SMESH_3D_Algo(hypId, studyId, gen)
  83 {
  84   MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
  85   _name = "Hexa_3D";
  86   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
  87   _requireShape = false;
  88   _compatibleHypothesis.push_back("ViscousLayers");
  89 }
  90
  91 //=============================================================================
  92 /*!
  93  * Destructor
  94  */
  95 //=============================================================================
  96
  97 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
  98 {
  99   MESSAGE("StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D");
 100 }
 101
 102 //=============================================================================
 103 /*!
 104  * Retrieves defined hypotheses
 105  */
 106 //=============================================================================
 107
 108 bool StdMeshers_Hexa_3D::CheckHypothesis
 109                          (SMESH_Mesh&                          aMesh,
 110                           const TopoDS_Shape&                  aShape,
 111                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 112 {
 113   // check nb of faces in the shape
 114 /*  PAL16229
 115   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
 116   int nbFaces = 0;
 117   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
 118     if ( ++nbFaces > 6 )
 119       break;
 120   if ( nbFaces != 6 )
 121     return false;
 122 */
 123
 124   _viscousLayersHyp = NULL;
 125
 126   const list<const SMESHDS_Hypothesis*>& hyps =
 127     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
 128   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
 129   if ( h == hyps.end())
 130   {
 131     aStatus = SMESH_Hypothesis::HYP_OK;
 132     return true;
 133   }
 134
 135   aStatus = HYP_OK;
 136   for ( ; h != hyps.end(); ++h )
 137   {
 138     string hypName = (*h)->GetName();
 139     if ( find( _compatibleHypothesis.begin(),_compatibleHypothesis.end(),hypName )
 140          != _compatibleHypothesis.end() )
 141     {
 142       _viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h );
 143     }
 144     else
 145     {
 146       aStatus = HYP_INCOMPATIBLE;
 147     }
 148   }
 149
 150   if ( !_viscousLayersHyp )
 151     aStatus = HYP_INCOMPATIBLE;
 152
 153   return aStatus == HYP_OK;
 154 }
 155
 156 namespace
 157 {
 158   //=============================================================================
 159
 160   // symbolic names of box sides
 161   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
 162
 163   // indices of FACE's of box sides in terms of SMESH_Block::TShapeID enum
 164   enum ESideIndex{ I_BOTTOM    = SMESH_Block::ID_Fxy0,
 165                    I_RIGHT     = SMESH_Block::ID_F1yz,
 166                    I_TOP       = SMESH_Block::ID_Fxy1,
 167                    I_LEFT      = SMESH_Block::ID_F0yz,
 168                    I_FRONT     = SMESH_Block::ID_Fx0z,
 169                    I_BACK      = SMESH_Block::ID_Fx1z,
 170                    I_UNDEFINED = SMESH_Block::ID_NONE
 171   };
 172   //=============================================================================
 173   /*!
 174    * \brief Container of nodes of structured mesh on a qudrangular geom FACE
 175    */
 176   struct _FaceGrid
 177   {
 178     // map of (node parameter on EDGE) to (column (vector) of nodes)
 179     TParam2ColumnMap _u2nodesMap;
 180
 181     // node column's taken form _u2nodesMap taking into account sub-shape orientation
 182     vector<TNodeColumn> _columns;
 183
 184     // geometry of a cube side
 185     TopoDS_Face _sideF;
 186     TopoDS_Edge _baseE;
 187
 188     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
 189     {
 190       return _columns[iCol][iRow];
 191     }
 192     gp_XYZ GetXYZ(int iCol, int iRow) const
 193     {
 194       return SMESH_MeshEditor::TNodeXYZ( GetNode( iCol, iRow ));
 195     }
 196   };
 197
 198   //================================================================================
 199   /*!
 200    * \brief Convertor of a pair of integers to a sole index
 201    */
 202   struct _Indexer
 203   {
 204     int _xSize, _ySize;
 205     _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
 206     int size() const { return _xSize * _ySize; }
 207     int operator()(const int x, const int y) const { return y * _xSize + x; }
 208   };
 209 }
 210
 211 //=============================================================================
 212 /*!
 213  * Generates hexahedron mesh on hexaedron like form using algorithm from
 214  * "Application de l'interpolation transfinie à la création de maillages
 215  *  C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
 216  *  et hexaedres déformés."
 217  * Alain PERONNET - 8 janvier 1999
 218  */
 219 //=============================================================================
 220
 221 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
 222                                  const TopoDS_Shape & aShape)// throw(SALOME_Exception)
 223 {
 224   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
 225   //Unexpect aCatch(SalomeException);
 226   MESSAGE("StdMeshers_Hexa_3D::Compute");
 227   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 228
 229   // 1) Shape verification
 230   // ----------------------
 231
 232   // shape must be a solid (or a shell) with 6 faces
 233   TopoDS_Shell shell;
 234   TopExp_Explorer exp(aShape,TopAbs_SHELL);
 235   if ( !exp.More() )
 236     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
 237   shell = TopoDS::Shell( exp.Current());
 238   if ( exp.Next(), exp.More() )
 239     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
 240
 241   TopTools_IndexedMapOfShape FF;
 242   TopExp::MapShapes( shell, TopAbs_FACE, FF);
 243   if ( FF.Extent() != 6)
 244   {
 245     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen);
 246     if ( !compositeHexa.Compute( aMesh, aShape ))
 247       return error( compositeHexa.GetComputeError() );
 248     return true;
 249   }
 250
 251   // Find sub-shapes of a cube
 252   TopTools_IndexedMapOfOrientedShape cubeSubShapes;
 253   TopoDS_Vertex V;
 254   if ( !SMESH_Block::FindBlockShapes( shell, V,V, cubeSubShapes ))
 255     return error(COMPERR_BAD_SHAPE, "Not a 6 sides cube");
 256
 257   // 2) make viscous layers
 258
 259   SMESH_ProxyMesh::Ptr proxymesh;
 260   if ( _viscousLayersHyp )
 261   {
 262     proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
 263     if ( !proxymesh )
 264       return false;
 265   }
 266
 267   // 3) Check presence of regular grid mesh on FACEs of the cube
 268   // ------------------------------------------------------------
 269
 270   // indices of FACEs of cube sides within cubeSubShapes
 271   const int sideIndex[6] = { I_BOTTOM, I_RIGHT, I_TOP, I_LEFT, I_FRONT, I_BACK };
 272   // indices of base EDGEs of cube sides within cubeSubShapes, corresponding to sideIndex
 273   const int baseEdgeIndex[6] = {
 274     SMESH_Block::ID_Ex00, // bottom side
 275     SMESH_Block::ID_E1y0, // right side
 276     SMESH_Block::ID_Ex01, // top side
 277     SMESH_Block::ID_E0y0, // left side
 278     SMESH_Block::ID_Ex00, // front side
 279     SMESH_Block::ID_Ex10  // back side
 280   };
 281
 282   // Load mesh of cube sides
 283
 284   _FaceGrid aCubeSide[ 6 ] ;
 285
 286   // tool creating quadratic elements if needed
 287   SMESH_MesherHelper helper (aMesh);
 288   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
 289
 290   for ( int i = 0; i < 6; ++i )
 291   {
 292     aCubeSide[i]._sideF = TopoDS::Face( cubeSubShapes( sideIndex[i] ));
 293     aCubeSide[i]._baseE = TopoDS::Edge( cubeSubShapes( baseEdgeIndex[i] ));
 294
 295     // assure correctness of node positions on _baseE
 296     if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( aCubeSide[i]._baseE ))
 297     {
 298       bool ok;
 299       helper.SetSubShape( aCubeSide[i]._baseE );
 300       SMDS_ElemIteratorPtr eIt = smDS->GetElements();
 301       while ( eIt->more() )
 302       {
 303         const SMDS_MeshElement* e = eIt->next();
 304         helper.GetNodeU( aCubeSide[i]._baseE, e->GetNode(0), e->GetNode(1), &ok);
 305         helper.GetNodeU( aCubeSide[i]._baseE, e->GetNode(1), e->GetNode(0), &ok);
 306       }
 307     }
 308
 309     // load grid
 310     if ( !helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap,
 311                                   aCubeSide[i]._sideF,
 312                                   aCubeSide[i]._baseE, meshDS, proxymesh.get() ))
 313     {
 314       SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 315       return error( err );
 316     }
 317
 318     // check if there are triangles on aCubeSide[i]._sideF
 319     if ( aMesh.NbTriangles() > 0 )
 320     {
 321       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( aCubeSide[i]._sideF ))
 322       {
 323         bool isAllQuad = true;
 324         SMDS_ElemIteratorPtr fIt = smDS->GetElements();
 325         while ( fIt->more() && isAllQuad )
 326         {
 327           const SMDS_MeshElement* f = fIt->next();
 328           isAllQuad = ( f->NbCornerNodes() == 4 );
 329         }
 330         if ( !isAllQuad )
 331         {
 332           SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
 333           return error( err );
 334         }
 335       }
 336     }
 337   }
 338
 339   // Orient loaded grids of cube sides along axis of the unitary cube coord system
 340   for ( int i = 0; i < 6; ++i )
 341   {
 342     bool reverse = false;
 343     if ( helper.GetSubShapeOri( shell.Oriented( TopAbs_FORWARD ),
 344                                 aCubeSide[i]._sideF ) == TopAbs_REVERSED )
 345       reverse = !reverse;
 346
 347     if ( helper.GetSubShapeOri( aCubeSide[i]._sideF.Oriented( TopAbs_FORWARD ),
 348                                 aCubeSide[i]._baseE ) == TopAbs_REVERSED )
 349       reverse = !reverse;
 350
 351     if ( sideIndex[i] == I_BOTTOM ||
 352          sideIndex[i] == I_LEFT   ||
 353          sideIndex[i] == I_BACK )
 354       reverse = !reverse;
 355
 356     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
 357
 358     int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
 359     int* pi = reverse ? &iRev : &iFwd;
 360     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
 361     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
 362       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
 363
 364     aCubeSide[i]._u2nodesMap.clear();
 365   }
 366
 367   if ( proxymesh )
 368     for ( int i = 0; i < 6; ++i )
 369       for ( unsigned j = 0; j < aCubeSide[i]._columns.size(); ++j)
 370         for ( unsigned k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
 371         {
 372           const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
 373           n = proxymesh->GetProxyNode( n );
 374         }
 375
 376   // 4) Create internal nodes of the cube
 377   // -------------------------------------
 378
 379   helper.SetSubShape( aShape );
 380   helper.SetElementsOnShape(true);
 381
 382   // shortcuts to sides
 383   _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
 384   _FaceGrid* fRight  = & aCubeSide[ B_RIGHT  ];
 385   _FaceGrid* fTop    = & aCubeSide[ B_TOP    ];
 386   _FaceGrid* fLeft   = & aCubeSide[ B_LEFT   ];
 387   _FaceGrid* fFront  = & aCubeSide[ B_FRONT  ];
 388   _FaceGrid* fBack   = & aCubeSide[ B_BACK   ];
 389
 390   // cube size measured in nb of nodes
 391   int x, xSize = fBottom->_columns.size() , X = xSize - 1;
 392   int y, ySize = fLeft->_columns.size()   , Y = ySize - 1;
 393   int z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
 394
 395   // columns of internal nodes "rising" from nodes of fBottom
 396   _Indexer colIndex( xSize, ySize );
 397   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
 398
 399   // fill node columns by front and back box sides
 400   for ( x = 0; x < xSize; ++x ) {
 401     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
 402     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
 403     column0.resize( zSize );
 404     column1.resize( zSize );
 405     for ( z = 0; z < zSize; ++z ) {
 406       column0[ z ] = fFront->GetNode( x, z );
 407       column1[ z ] = fBack ->GetNode( x, z );
 408     }
 409   }
 410   // fill node columns by left and right box sides
 411   for ( y = 1; y < ySize-1; ++y ) {
 412     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
 413     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
 414     column0.resize( zSize );
 415     column1.resize( zSize );
 416     for ( z = 0; z < zSize; ++z ) {
 417       column0[ z ] = fLeft ->GetNode( y, z );
 418       column1[ z ] = fRight->GetNode( y, z );
 419     }
 420   }
 421   // get nodes from top and bottom box sides
 422   for ( x = 1; x < xSize-1; ++x ) {
 423     for ( y = 1; y < ySize-1; ++y ) {
 424       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 425       column.resize( zSize );
 426       column.front() = fBottom->GetNode( x, y );
 427       column.back()  = fTop   ->GetNode( x, y );
 428     }
 429   }
 430
 431   // projection points of the internal node on cube sub-shapes by which
 432   // coordinates of the internal node are computed
 433   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
 434
 435   // projections on vertices are constant
 436   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
 437   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
 438   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
 439   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
 440   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
 441   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
 442   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
 443   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
 444
 445   for ( x = 1; x < xSize-1; ++x )
 446   {
 447     gp_XYZ params; // normalized parameters of internal node within a unit box
 448     params.SetCoord( 1, x / double(X) );
 449     for ( y = 1; y < ySize-1; ++y )
 450     {
 451       params.SetCoord( 2, y / double(Y) );
 452       // a column to fill in during z loop
 453       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
 454       // projection points on horizontal edges
 455       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
 456       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
 457       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
 458       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
 459       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
 460       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
 461       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
 462       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
 463       // projection points on horizontal faces
 464       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
 465       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
 466       for ( z = 1; z < zSize-1; ++z ) // z loop
 467       {
 468         params.SetCoord( 3, z / double(Z) );
 469         // projection points on vertical edges
 470         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
 471         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
 472         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
 473         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
 474         // projection points on vertical faces
 475         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
 476         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
 477         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
 478         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
 479
 480         // compute internal node coordinates
 481         gp_XYZ coords;
 482         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
 483         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
 484
 485       }
 486     }
 487   }
 488
 489   // side data no more needed, free memory
 490   for ( int i = 0; i < 6; ++i )
 491     aCubeSide[i]._columns.clear();
 492
 493   // 5) Create hexahedrons
 494   // ---------------------
 495
 496   for ( x = 0; x < xSize-1; ++x ) {
 497     for ( y = 0; y < ySize-1; ++y ) {
 498       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
 499       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
 500       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
 501       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
 502       for ( z = 0; z < zSize-1; ++z )
 503       {
 504         // bottom face normal of a hexa mush point outside the volume
 505         helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
 506                          col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
 507       }
 508     }
 509   }
 510   return true;
 511 }
 512
 513 //=============================================================================
 514 /*!
 515  *  Evaluate
 516  */
 517 //=============================================================================
 518
 519 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
 520                                   const TopoDS_Shape & aShape,
 521                                   MapShapeNbElems& aResMap)
 522 {
 523   vector < SMESH_subMesh * >meshFaces;
 524   TopTools_SequenceOfShape aFaces;
 525   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
 526     aFaces.Append(exp.Current());
 527     SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 528     ASSERT(aSubMesh);
 529     meshFaces.push_back(aSubMesh);
 530   }
 531   if (meshFaces.size() != 6) {
 532     //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
 533     static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
 534     return compositeHexa.Evaluate(aMesh, aShape, aResMap);
 535   }
 536
 537   int i = 0;
 538   for(; i<6; i++) {
 539     //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
 540     TopoDS_Shape aFace = aFaces.Value(i+1);
 541     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
 542     if( !algo ) {
 543       std::vector<int> aResVec(SMDSEntity_Last);
 544       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 545       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 546       aResMap.insert(std::make_pair(sm,aResVec));
 547       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 548       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 549       return false;
 550     }
 551     string algoName = algo->GetName();
 552     bool isAllQuad = false;
 553     if (algoName == "Quadrangle_2D") {
 554       MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
 555       if( anIt == aResMap.end() ) continue;
 556       std::vector<int> aVec = (*anIt).second;
 557       int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 558       if( nbtri == 0 )
 559         isAllQuad = true;
 560     }
 561     if ( ! isAllQuad ) {
 562       return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
 563     }
 564   }
 565
 566   // find number of 1d elems for 1 face
 567   int nb1d = 0;
 568   TopTools_MapOfShape Edges1;
 569   bool IsQuadratic = false;
 570   bool IsFirst = true;
 571   for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
 572     Edges1.Add(exp.Current());
 573     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
 574     if( sm ) {
 575       MapShapeNbElemsItr anIt = aResMap.find(sm);
 576       if( anIt == aResMap.end() ) continue;
 577       std::vector<int> aVec = (*anIt).second;
 578       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 579       if(IsFirst) {
 580         IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
 581         IsFirst = false;
 582       }
 583     }
 584   }
 585   // find face opposite to 1 face
 586   int OppNum = 0;
 587   for(i=2; i<=6; i++) {
 588     bool IsOpposite = true;
 589     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
 590       if( Edges1.Contains(exp.Current()) ) {
 591         IsOpposite = false;
 592         break;
 593       }
 594     }
 595     if(IsOpposite) {
 596       OppNum = i;
 597       break;
 598     }
 599   }
 600   // find number of 2d elems on side faces
 601   int nb2d = 0;
 602   for(i=2; i<=6; i++) {
 603     if( i == OppNum ) continue;
 604     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
 605     if( anIt == aResMap.end() ) continue;
 606     std::vector<int> aVec = (*anIt).second;
 607     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 608   }
 609
 610   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
 611   std::vector<int> aVec = (*anIt).second;
 612   int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 613   int nb0d_face0 = aVec[SMDSEntity_Node];
 614
 615   std::vector<int> aResVec(SMDSEntity_Last);
 616   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 617   if(IsQuadratic) {
 618     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
 619     int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
 620     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
 621   }
 622   else {
 623     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
 624     aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
 625   }
 626   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
 627   aResMap.insert(std::make_pair(sm,aResVec));
 628
 629   return true;
 630 }
 631
 632 //================================================================================
 633 /*!
 634  * \brief Computes hexahedral mesh from 2D mesh of block
 635  */
 636 //================================================================================
 637
 638 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
 639 {
 640   static StdMeshers_HexaFromSkin_3D * algo = 0;
 641   if ( !algo ) {
 642     SMESH_Gen* gen = aMesh.GetGen();
 643     algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), 0, gen );
 644   }
 645   algo->InitComputeError();
 646   algo->Compute( aMesh, aHelper );
 647   return error( algo->GetComputeError());
 648 }
 649
 650 //=======================================================================
 651 //function : ComputePentahedralMesh
 652 //purpose  :
 653 //=======================================================================
 654
 655 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &          aMesh,
 656                                              const TopoDS_Shape &  aShape,
 657                                              SMESH_ProxyMesh*      proxyMesh)
 658 {
 659   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
 660   if ( proxyMesh )
 661   {
 662     err->myName = COMPERR_BAD_INPUT_MESH;
 663     err->myComment = "Can't build pentahedral mesh on viscous layers";
 664     return err;
 665   }
 666   bool bOK;
 667   StdMeshers_Penta_3D anAlgo;
 668   //
 669   bOK=anAlgo.Compute(aMesh, aShape);
 670   //
 671   err = anAlgo.GetComputeError();
 672   //
 673   if ( !bOK && anAlgo.ErrorStatus() == 5 )
 674   {
 675     static StdMeshers_Prism_3D * aPrism3D = 0;
 676     if ( !aPrism3D ) {
 677       SMESH_Gen* gen = aMesh.GetGen();
 678       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
 679     }
 680     SMESH_Hypothesis::Hypothesis_Status aStatus;
 681     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
 682       aPrism3D->InitComputeError();
 683       bOK = aPrism3D->Compute( aMesh, aShape );
 684       err = aPrism3D->GetComputeError();
 685     }
 686   }
 687   return err;
 688 }
 689
 690
 691 //=======================================================================
 692 //function : EvaluatePentahedralMesh
 693 //purpose  :
 694 //=======================================================================
 695
 696 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
 697                              const TopoDS_Shape & aShape,
 698                              MapShapeNbElems& aResMap)
 699 {
 700   StdMeshers_Penta_3D anAlgo;
 701   bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
 702
 703   //err = anAlgo.GetComputeError();
 704   //if ( !bOK && anAlgo.ErrorStatus() == 5 )
 705   if( !bOK ) {
 706     static StdMeshers_Prism_3D * aPrism3D = 0;
 707     if ( !aPrism3D ) {
 708       SMESH_Gen* gen = aMesh.GetGen();
 709       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
 710     }
 711     SMESH_Hypothesis::Hypothesis_Status aStatus;
 712     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
 713       return aPrism3D->Evaluate(aMesh, aShape, aResMap);
 714     }
 715   }
 716
 717   return bOK;
 718 }