src/StdMeshers/StdMeshers_Hexa_3D.cxx

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