src/StdMeshers/StdMeshers_Hexa_3D.cxx

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