src/StdMeshers/StdMeshers_Hexa_3D.cxx

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