src/StdMeshers/StdMeshers_Hexa_3D.cxx

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