Salome HOME
Fix compilation errors under Windows: VS-2017
[modules/smesh.git] / src / StdMeshers / StdMeshers_Hexa_3D.cxx
1 // Copyright (C) 2007-2020  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 "SMDS_MeshNode.hxx"
32 #include "SMESH_Comment.hxx"
33 #include "SMESH_Gen.hxx"
34 #include "SMESH_Mesh.hxx"
35 #include "SMESH_MesherHelper.hxx"
36 #include "SMESH_subMesh.hxx"
37 #include "StdMeshers_BlockRenumber.hxx"
38 #include "StdMeshers_CompositeHexa_3D.hxx"
39 #include "StdMeshers_FaceSide.hxx"
40 #include "StdMeshers_HexaFromSkin_3D.hxx"
41 #include "StdMeshers_Penta_3D.hxx"
42 #include "StdMeshers_Prism_3D.hxx"
43 #include "StdMeshers_Quadrangle_2D.hxx"
44 #include "StdMeshers_ViscousLayers.hxx"
45
46 #include <BRep_Tool.hxx>
47 #include <Bnd_B3d.hxx>
48 #include <TopExp.hxx>
49 #include <TopExp_Explorer.hxx>
50 #include <TopTools_MapOfShape.hxx>
51 #include <TopTools_SequenceOfShape.hxx>
52 #include <TopoDS.hxx>
53
54 #include "utilities.h"
55 #include "Utils_ExceptHandlers.hxx"
56
57 #include <cstddef>
58
59 #include <numeric>
60
61 typedef SMESH_Comment TComm;
62
63 using namespace std;
64
65 static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
66                                                     const TopoDS_Shape &,
67                                                     SMESH_ProxyMesh* proxyMesh=0);
68
69 static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
70                                     MapShapeNbElems &);
71
72 //=============================================================================
73 /*!
74  * Constructor
75  */
76 //=============================================================================
77
78 StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, SMESH_Gen * gen)
79   :SMESH_3D_Algo(hypId, gen)
80 {
81   _name = "Hexa_3D";
82   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
83   _requireShape = false;
84   _compatibleHypothesis.push_back("ViscousLayers");
85   _compatibleHypothesis.push_back("BlockRenumber");
86   _quadAlgo = new StdMeshers_Quadrangle_2D( gen->GetANewId(), _gen );
87 }
88
89 //=============================================================================
90 /*!
91  * Destructor
92  */
93 //=============================================================================
94
95 StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
96 {
97   delete _quadAlgo;
98   _quadAlgo = 0;
99 }
100
101 //=============================================================================
102 /*!
103  * Retrieves defined hypotheses
104  */
105 //=============================================================================
106
107 bool StdMeshers_Hexa_3D::CheckHypothesis
108                          (SMESH_Mesh&                          aMesh,
109                           const TopoDS_Shape&                  aShape,
110                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
111 {
112   // check nb of faces in the shape
113 /*  PAL16229
114   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
115   int nbFaces = 0;
116   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
117     if ( ++nbFaces > 6 )
118       break;
119   if ( nbFaces != 6 )
120     return false;
121 */
122
123   _viscousLayersHyp = nullptr;
124   _blockRenumberHyp = nullptr;
125
126   const list<const SMESHDS_Hypothesis*>& hyps =
127     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliary=*/false);
128   list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
129   if ( h == hyps.end())
130   {
131     aStatus = SMESH_Hypothesis::HYP_OK;
132     return true;
133   }
134
135   // only StdMeshers_ViscousLayers can be used
136   aStatus = HYP_OK;
137   for ( ; h != hyps.end(); ++h )
138   {
139     if ( !_viscousLayersHyp &&
140          (_viscousLayersHyp = dynamic_cast< const StdMeshers_ViscousLayers*> ( *h )))
141       continue;
142     if ( !_blockRenumberHyp &&
143          (_blockRenumberHyp = dynamic_cast< const StdMeshers_BlockRenumber*> ( *h )))
144       continue;
145     break;
146   }
147   if ((int) hyps.size() != (bool)_viscousLayersHyp + (bool)_blockRenumberHyp )
148     aStatus = HYP_INCOMPATIBLE;
149   else
150   {
151     if ( _viscousLayersHyp )
152       if ( !error( _viscousLayersHyp->CheckHypothesis( aMesh, aShape, aStatus )))
153         aStatus = HYP_BAD_PARAMETER;
154
155     if ( _blockRenumberHyp && aStatus == HYP_OK )
156       error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape ));
157   }
158
159   return aStatus == HYP_OK;
160 }
161
162 namespace
163 {
164   //=============================================================================
165
166   typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
167   typedef std::vector<gp_XYZ>                 TXYZColumn;
168
169   // symbolic names of box sides
170   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
171
172   // symbolic names of sides of quadrangle
173   enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
174
175   enum EAxes{ COO_X=1, COO_Y, COO_Z };
176
177   //=============================================================================
178   /*!
179    * \brief Container of nodes of structured mesh on a qudrangular geom FACE
180    */
181   struct _FaceGrid
182   {
183     // face sides
184     FaceQuadStructPtr _quad;
185
186     // map of (node parameter on EDGE) to (column (vector) of nodes)
187     TParam2ColumnMap _u2nodesMap;
188
189     // node column's taken from _u2nodesMap taking into account sub-shape orientation
190     vector<TNodeColumn> _columns;
191
192     // columns of normalized parameters of nodes within the unitary cube
193     vector<TXYZColumn> _ijkColumns;
194
195     // geometry of a cube side
196     TopoDS_Face _sideF;
197
198     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
199     {
200       return _columns[iCol][iRow];
201     }
202     gp_XYZ GetXYZ(int iCol, int iRow) const
203     {
204       return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
205     }
206     gp_XYZ& GetIJK(int iCol, int iRow)
207     {
208       return _ijkColumns[iCol][iRow];
209     }
210   };
211
212   //================================================================================
213   /*!
214    * \brief Converter of a pair of integers to a sole index
215    */
216   struct _Indexer
217   {
218     int _xSize, _ySize;
219     _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
220     int size() const { return _xSize * _ySize; }
221     int operator()(const int x, const int y) const { return y * _xSize + x; }
222   };
223
224   //================================================================================
225   /*!
226    * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
227    *  the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
228    */
229   FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSidePtr side,
230                                        FaceQuadStructPtr      quad[ 6 ])
231   {
232     FaceQuadStructPtr foundQuad;
233     for ( int i = 1; i < 6; ++i )
234     {
235       if ( !quad[i] ) continue;
236       for ( size_t iS = 0; iS < quad[i]->side.size(); ++iS )
237       {
238         const StdMeshers_FaceSidePtr side2 = quad[i]->side[iS];
239         if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
240               side->FirstVertex().IsSame( side2->LastVertex() ))
241             &&
242             ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
243               side->LastVertex().IsSame( side2->LastVertex() ))
244             )
245         {
246           if ( iS != Q_BOTTOM )
247           {
248             vector< FaceQuadStruct::Side > newSides;
249             for ( size_t j = iS; j < quad[i]->side.size(); ++j )
250               newSides.push_back( quad[i]->side[j] );
251             for ( size_t j = 0; j < iS; ++j )
252               newSides.push_back( quad[i]->side[j] );
253             quad[i]->side.swap( newSides );
254           }
255           foundQuad.swap(quad[i]);
256           return foundQuad;
257         }
258       }
259     }
260     return foundQuad;
261   }
262
263   //================================================================================
264   /*!
265    * \brief Put quads to aCubeSide in the order of enum EBoxSides
266    */
267   //================================================================================
268
269   bool arrangeQuads( FaceQuadStructPtr quad[ 6 ], _FaceGrid aCubeSide[ 6 ], bool reverseBottom )
270   {
271     swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
272     if ( reverseBottom )
273       swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the bottom normal inside cube
274             aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
275
276     aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
277     aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
278     aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP   ], quad );
279     aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT  ], quad );
280     if ( aCubeSide[B_FRONT ]._quad )
281       aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
282
283     for ( int i = 1; i < 6; ++i )
284       if ( !aCubeSide[i]._quad )
285         return false;
286     return true;
287   }
288
289   //================================================================================
290   /*!
291    * \brief Rearrange block sides according to StdMeshers_BlockRenumber hypothesis
292    */
293   //================================================================================
294
295   bool arrangeForRenumber( _FaceGrid      blockSide[ 6 ],
296                            TopoDS_Vertex& v000,
297                            TopoDS_Vertex& v001 )
298   {
299     std::swap( blockSide[B_BOTTOM]._quad->side[ Q_RIGHT],// restore after arrangeQuads()
300                blockSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
301
302     // find v000
303     TopTools_MapOfShape cornerVertices;
304     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_BOTTOM].grid->LastVertex()  );
305     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_BOTTOM].grid->FirstVertex() );
306     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_TOP   ].grid->LastVertex()  );
307     cornerVertices.Add(  blockSide[B_BOTTOM]._quad->side[Q_TOP   ].grid->FirstVertex() );
308     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_BOTTOM].grid->FirstVertex() );
309     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_BOTTOM].grid->LastVertex()  );
310     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_TOP   ].grid->FirstVertex() );
311     cornerVertices.Add(  blockSide[B_TOP   ]._quad->side[Q_TOP   ].grid->LastVertex()  );
312
313     if ( v000.IsNull() )
314     {
315       // block CS is not defined;
316       // renumber only if the block has an edge parallel to an axis of global CS
317
318       v000 = StdMeshers_RenumberHelper::GetVertex000( cornerVertices );
319     }
320
321     Bnd_B3d bbox;
322     for ( auto it = cornerVertices.cbegin(); it != cornerVertices.cend(); ++it )
323       bbox.Add( BRep_Tool::Pnt( TopoDS::Vertex( *it )));
324     double tol = 1e-5 * Sqrt( bbox.SquareExtent() );
325
326     // get block edges starting at v000
327
328     std::vector< StdMeshers_FaceSidePtr > edgesAtV000;
329     std::vector< gp_Vec >                 edgeDir;
330     std::vector< int >                    iParallel; // 0 - none, 1 - X, 2 - Y, 3 - Z
331     TopTools_MapOfShape                   lastVertices;
332     for ( int iQ = 0; iQ < 6; ++iQ )
333     {
334       FaceQuadStructPtr quad = blockSide[iQ]._quad;
335       for ( size_t iS = 0; iS < quad->side.size() &&  edgesAtV000.size() < 3; ++iS )
336       {
337         StdMeshers_FaceSidePtr edge = quad->side[iS];
338         TopoDS_Vertex v1 = edge->FirstVertex(), v2 = edge->LastVertex();
339         if (( v1.IsSame( v000 ) && !lastVertices.Contains( v2 )) ||
340             ( v2.IsSame( v000 ) && !lastVertices.Contains( v1 )))
341         {
342           bool reverse = v2.IsSame( v000 );
343           if ( reverse )
344             std::swap( v1, v2 );
345           lastVertices.Add( v2 );
346
347           edgesAtV000.push_back( edge );
348
349           gp_Pnt pf = BRep_Tool::Pnt( v1 );
350           gp_Pnt pl = BRep_Tool::Pnt( v2 );
351           gp_Vec vec( pf, pl );
352           edgeDir.push_back( vec );
353
354           iParallel.push_back( 0 );
355           if ( !v001.IsNull() )
356           {
357             if ( v001.IsSame( v2 ))
358               iParallel.back() = 3;
359           }
360           else
361           {
362             bool isStraight = true;
363             for ( int iE = 0; iE < edge->NbEdges() &&  isStraight; ++iE )
364               isStraight = SMESH_Algo::IsStraight( edge->Edge( iE ));
365
366             // is parallel to a GCS axis?
367             if ( isStraight )
368             {
369               int nbDiff = (( Abs( vec.X() ) > tol ) +
370                             ( Abs( vec.Y() ) > tol ) +
371                             ( Abs( vec.Z() ) > tol ) );
372               if ( nbDiff == 1 )
373                 iParallel.back() = ( Abs( vec.X() ) > tol ) ? 1 : ( Abs( vec.Y() ) > tol ) ? 2 : 3;
374             }
375             else
376             {
377               edgeDir.back() = gp_Vec( pf, edge->Value3d( reverse ? 0.99 : 0.01 ));
378             }
379           }
380         }
381       }
382     }
383     if ( std::accumulate( iParallel.begin(), iParallel.end(), 0 ) == 0 )
384       return false;
385
386     // find edge OZ and edge OX
387     StdMeshers_FaceSidePtr edgeOZ, edgeOX;
388     auto iZIt = std::find( iParallel.begin(), iParallel.end(), 3 );
389     if ( iZIt != iParallel.end() )
390     {
391       int i = std::distance( iParallel.begin(), iZIt );
392       edgeOZ = edgesAtV000[ i ];
393       int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() );
394       int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() );
395       if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 )
396         std::swap( iE1, iE2 );
397       edgeOX = edgesAtV000[ iE1 ];
398     }
399     else
400     {
401       for ( size_t i = 0; i < edgesAtV000.size(); ++i )
402       {
403         if ( !iParallel[ i ] )
404           continue;
405         int iE1 = SMESH_MesherHelper::WrapIndex( i + 1, edgesAtV000.size() );
406         int iE2 = SMESH_MesherHelper::WrapIndex( i + 2, edgesAtV000.size() );
407         if (( edgeDir[ iE1 ] ^ edgeDir[ iE2 ] ) * edgeDir[ i ] < 0 )
408           std::swap( iE1, iE2 );
409         edgeOZ = edgesAtV000[ iParallel[i] == 1 ? iE2 : iE1 ];
410         edgeOX = edgesAtV000[ iParallel[i] == 1 ? i : iE1 ];
411         break;
412       }
413     }
414
415     if ( !edgeOZ || !edgeOX )
416       return false;
417
418     TopoDS_Vertex v100 = edgeOX->LastVertex();
419     if ( v100.IsSame( v000 ))
420       v100 = edgeOX->FirstVertex();
421
422     // Find the left quad, one including v000 but not v100
423
424     for ( int iQ = 0; iQ < 6; ++iQ )
425     {
426       FaceQuadStructPtr quad = blockSide[iQ]._quad;
427       bool hasV000 = false, hasV100 = false;
428       for ( size_t iS = 0; iS < quad->side.size(); ++iS )
429       {
430         StdMeshers_FaceSidePtr edge = quad->side[iS];
431         if ( edge->FirstVertex().IsSame( v000 ) || edge->LastVertex().IsSame( v000 ))
432           hasV000 = true;
433         if ( edge->FirstVertex().IsSame( v100 ) || edge->LastVertex().IsSame( v100 ))
434           hasV100 = true;
435       }
436       if ( hasV000 && !hasV100 )
437       {
438         // orient the left quad
439         for ( int i = 0; i < 4; ++i )
440         {
441           if ( quad->side[Q_BOTTOM].grid->Edge(0).IsSame( edgeOZ->Edge(0) ))
442             break;
443           quad->shift( 1, true );
444         }
445
446         FaceQuadStructPtr quads[ 6 ];
447         quads[0].swap( blockSide[iQ]._quad );
448         for ( int i = 1, j = 0; i < 6; ++i, ++j )
449           if ( blockSide[ j ]._quad )
450             quads[ i ].swap( blockSide[ j ]._quad );
451           else
452             --i;
453
454         return arrangeQuads( quads, blockSide, false/* true*/ );
455       }
456     }
457     return false;
458   }
459
460   //================================================================================
461   /*!
462    * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
463    */
464   //================================================================================
465
466   bool beginsAtSide( const _FaceGrid&     sideGrid1,
467                      const _FaceGrid&     sideGrid2,
468                      SMESH_ProxyMesh::Ptr proxymesh )
469   {
470     const TNodeColumn& col0  = sideGrid2._u2nodesMap.begin()->second;
471     const TNodeColumn& col1  = sideGrid2._u2nodesMap.rbegin()->second;
472     const SMDS_MeshNode* n00 = col0.front();
473     const SMDS_MeshNode* n01 = col0.back();
474     const SMDS_MeshNode* n10 = col1.front();
475     const SMDS_MeshNode* n11 = col1.back();
476     const SMDS_MeshNode* n = (sideGrid1._u2nodesMap.begin()->second)[0];
477     if ( proxymesh )
478     {
479       n00 = proxymesh->GetProxyNode( n00 );
480       n10 = proxymesh->GetProxyNode( n10 );
481       n01 = proxymesh->GetProxyNode( n01 );
482       n11 = proxymesh->GetProxyNode( n11 );
483       n   = proxymesh->GetProxyNode( n );
484     }
485     return ( n == n00 || n == n01 || n == n10 || n == n11 );
486   }
487
488   //================================================================================
489   /*!
490    * \brief Fill in _FaceGrid::_ijkColumns
491    *  \param [in,out] fg - a _FaceGrid
492    *  \param [in] i1 - coordinate index along _columns
493    *  \param [in] i2 - coordinate index along _columns[i]
494    *  \param [in] v3 - value of the constant parameter
495    */
496   //================================================================================
497
498   void computeIJK( _FaceGrid& fg, int i1, int i2, double v3 )
499   {
500     gp_XYZ ijk( v3, v3, v3 );
501     const size_t nbCol = fg._columns.size();
502     const size_t nbRow = fg._columns[0].size();
503
504     fg._ijkColumns.resize( nbCol );
505     for ( size_t i = 0; i < nbCol; ++i )
506       fg._ijkColumns[ i ].resize( nbRow, ijk );
507
508     vector< double > len( nbRow );
509     len[0] = 0;
510     for ( size_t i = 0; i < nbCol; ++i )
511     {
512       gp_Pnt pPrev = fg.GetXYZ( i, 0 );
513       for ( size_t j = 1; j < nbRow; ++j )
514       {
515         gp_Pnt p = fg.GetXYZ( i, j );
516         len[ j ] = len[ j-1 ] + p.Distance( pPrev );
517         pPrev = p;
518       }
519       for ( size_t j = 0; j < nbRow; ++j )
520         fg.GetIJK( i, j ).SetCoord( i2, len[ j ]/len.back() );
521     }
522
523     len.resize( nbCol );
524     for ( size_t j = 0; j < nbRow; ++j )
525     {
526       gp_Pnt pPrev = fg.GetXYZ( 0, j );
527       for ( size_t i = 1; i < nbCol; ++i )
528       {
529         gp_Pnt p = fg.GetXYZ( i, j );
530         len[ i ] = len[ i-1 ] + p.Distance( pPrev );
531         pPrev = p;
532       }
533       for ( size_t i = 0; i < nbCol; ++i )
534         fg.GetIJK( i, j ).SetCoord( i1, len[ i ]/len.back() );
535     }
536   }
537 }
538
539 //=============================================================================
540 /*!
541  * Generates hexahedron mesh on hexaedron like form using algorithm from
542  * "Application de l'interpolation transfinie � la cr�ation de maillages
543  *  C0 ou G1 continus sur des triangles, quadrangles, tetraedres, pentaedres
544  *  et hexaedres d�form�s."
545  * Alain PERONNET - 8 janvier 1999
546  */
547 //=============================================================================
548
549 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
550                                  const TopoDS_Shape & aShape)
551 {
552   // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
553   //Unexpect aCatch(SalomeException);
554   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
555
556   // Shape verification
557   // ----------------------
558
559   // shape must be a solid (or a shell) with 6 faces
560   TopExp_Explorer exp(aShape,TopAbs_SHELL);
561   if ( !exp.More() )
562     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
563   if ( exp.Next(), exp.More() )
564     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
565
566   TopTools_IndexedMapOfShape FF, EE;
567   TopExp::MapShapes( aShape, TopAbs_FACE, FF);
568   if ( FF.Extent() != 6)
569   {
570     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), _gen);
571     compositeHexa.SetHypothesis( _blockRenumberHyp );
572     if ( !compositeHexa.Compute( aMesh, aShape ))
573       return error( compositeHexa.GetComputeError() );
574
575     return _blockRenumberHyp ? error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape )) : true;
576   }
577
578   // Find sides of a cube
579   // ---------------------
580
581   // tool creating quadratic elements if needed
582   SMESH_MesherHelper helper (aMesh);
583   _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
584
585   TopExp::MapShapes( aShape, TopAbs_EDGE, EE );
586   SMESH_MesherHelper* faceHelper = ( EE.Size() == 12 ) ? 0 : &helper;
587
588   FaceQuadStructPtr quad[ 6 ];
589   for ( int i = 0; i < 6; ++i )
590   {
591     if ( faceHelper )
592       faceHelper->SetSubShape( FF( i+1 ));
593     if ( !( quad[i] = FaceQuadStructPtr( _quadAlgo->CheckNbEdges( aMesh, FF( i+1 ),
594                                                                   /*considerMesh=*/true,
595                                                                   faceHelper))))
596       return error( _quadAlgo->GetComputeError() );
597     if ( quad[i]->side.size() != 4 )
598       return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
599   }
600
601   // put quads in a proper order
602   _FaceGrid aCubeSide[ 6 ];
603   if ( !arrangeQuads( quad, aCubeSide, true ))
604     return error( COMPERR_BAD_SHAPE );
605
606
607   // Make viscous layers
608   // --------------------
609
610   SMESH_ProxyMesh::Ptr proxymesh;
611   if ( _viscousLayersHyp )
612   {
613     proxymesh = _viscousLayersHyp->Compute( aMesh, aShape, /*makeN2NMap=*/ true );
614     if ( !proxymesh )
615       return false;
616   }
617
618   // Check if there are triangles on cube sides
619   // -------------------------------------------
620
621   if ( aMesh.NbTriangles() > 0 )
622   {
623     for ( int i = 0; i < 6; ++i )
624     {
625       const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
626       const SMESHDS_SubMesh* smDS =
627         proxymesh ? proxymesh->GetSubMesh( sideF ) : meshDS->MeshElements( sideF );
628       if ( !SMESH_MesherHelper::IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE,
629                                                     /*nullSubMeshRes=*/false ))
630       {
631         SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
632         return error( err );
633       }
634     }
635   }
636
637   // Arrange sides according to _blockRenumberHyp
638   bool toRenumber = _blockRenumberHyp;
639   if ( toRenumber )
640   {
641     TopoDS_Vertex v000, v001;
642     _blockRenumberHyp->IsSolidIncluded( aMesh, aShape, v000, v001 );
643
644     toRenumber = arrangeForRenumber( aCubeSide, v000, v001 );
645
646     if ( toRenumber )
647     {
648       meshDS->Modified();
649       meshDS->CompactMesh(); // remove numbering holes
650     }
651   }
652
653   // Check presence of regular grid mesh on FACEs of the cube
654   // ------------------------------------------------------------
655
656   for ( int i = 0; i < 6; ++i )
657   {
658     const TopoDS_Face& F = aCubeSide[i]._quad->face;
659     StdMeshers_FaceSidePtr baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
660     list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
661
662     // assure correctness of node positions on baseE:
663     // helper.GetNodeU() will fix positions if they are wrong
664     helper.ToFixNodeParameters( true );
665     for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
666     {
667       const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
668       if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
669       {
670         bool ok;
671         helper.SetSubShape( baseE );
672         SMDS_ElemIteratorPtr eIt = smDS->GetElements();
673         while ( eIt->more() )
674         {
675           const SMDS_MeshElement* e = eIt->next();
676           // expect problems on a composite side
677           try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
678           catch (...) {}
679           try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
680           catch (...) {}
681         }
682       }
683     }
684
685     // load grid
686     bool ok =
687       helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
688     if ( ok )
689     {
690       // check if the loaded grid corresponds to nb of quadrangles on the FACE
691       const SMESHDS_SubMesh* faceSubMesh =
692         proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
693       const int nbQuads = faceSubMesh->NbElements();
694       const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
695       const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
696       ok = ( nbQuads == nbHor * nbVer );
697     }
698     if ( !ok )
699     {
700       SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
701       return error( err );
702     }
703   }
704
705   // Orient loaded grids of cube sides along axis of the unitary cube coord system
706   bool isReverse[6];
707   isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ], proxymesh );
708   isReverse[B_TOP   ] = beginsAtSide( aCubeSide[B_TOP   ], aCubeSide[B_RIGHT ], proxymesh );
709   isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ], proxymesh );
710   isReverse[B_BACK  ] = beginsAtSide( aCubeSide[B_BACK  ], aCubeSide[B_RIGHT ], proxymesh );
711   isReverse[B_LEFT  ] = beginsAtSide( aCubeSide[B_LEFT  ], aCubeSide[B_BACK  ], proxymesh );
712   isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK  ], proxymesh );
713   for ( int i = 0; i < 6; ++i )
714   {
715     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
716
717     size_t iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
718     size_t*  pi = isReverse[i] ? &iRev : &iFwd;
719     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
720     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
721       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );
722
723     aCubeSide[i]._u2nodesMap.clear();
724   }
725
726   if ( proxymesh )
727     for ( int i = 0; i < 6; ++i )
728       for ( size_t j = 0; j < aCubeSide[i]._columns.size(); ++j)
729         for ( size_t k = 0; k < aCubeSide[i]._columns[j].size(); ++k)
730         {
731           const SMDS_MeshNode* & n = aCubeSide[i]._columns[j][k];
732           n = proxymesh->GetProxyNode( n );
733         }
734
735   // 4) Create internal nodes of the cube
736   // -------------------------------------
737
738   helper.SetSubShape( aShape );
739   helper.SetElementsOnShape(true);
740
741   // shortcuts to sides
742   _FaceGrid* fBottom = & aCubeSide[ B_BOTTOM ];
743   _FaceGrid* fRight  = & aCubeSide[ B_RIGHT  ];
744   _FaceGrid* fTop    = & aCubeSide[ B_TOP    ];
745   _FaceGrid* fLeft   = & aCubeSide[ B_LEFT   ];
746   _FaceGrid* fFront  = & aCubeSide[ B_FRONT  ];
747   _FaceGrid* fBack   = & aCubeSide[ B_BACK   ];
748
749   // cube size measured in nb of nodes
750   size_t x, xSize = fBottom->_columns.size() , X = xSize - 1;
751   size_t y, ySize = fLeft->_columns.size()   , Y = ySize - 1;
752   size_t z, zSize = fLeft->_columns[0].size(), Z = zSize - 1;
753
754   // check sharing of FACEs (IPAL54417)
755   if ( fFront ->_columns.size()    != xSize ||
756        fBack  ->_columns.size()    != xSize ||
757        fTop   ->_columns.size()    != xSize ||
758
759        fRight ->_columns.size()    != ySize ||
760        fTop   ->_columns[0].size() != ySize ||
761        fBottom->_columns[0].size() != ySize ||
762
763        fRight ->_columns[0].size() != zSize ||
764        fFront ->_columns[0].size() != zSize ||
765        fBack  ->_columns[0].size() != zSize )
766     return error( COMPERR_BAD_SHAPE, "Not sewed faces" );
767
768   // columns of internal nodes "rising" from nodes of fBottom
769   _Indexer colIndex( xSize, ySize );
770   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
771
772   // fill node columns by front and back box sides
773   for ( x = 0; x < xSize; ++x ) {
774     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
775     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
776     column0.resize( zSize );
777     column1.resize( zSize );
778     for ( z = 0; z < zSize; ++z ) {
779       column0[ z ] = fFront->GetNode( x, z );
780       column1[ z ] = fBack ->GetNode( x, z );
781     }
782   }
783   // fill node columns by left and right box sides
784   for ( y = 1; y < ySize-1; ++y ) {
785     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
786     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
787     column0.resize( zSize );
788     column1.resize( zSize );
789     for ( z = 0; z < zSize; ++z ) {
790       column0[ z ] = fLeft ->GetNode( y, z );
791       column1[ z ] = fRight->GetNode( y, z );
792     }
793   }
794   // get nodes from top and bottom box sides
795   for ( x = 1; x < xSize-1; ++x ) {
796     for ( y = 1; y < ySize-1; ++y ) {
797       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
798       column.resize( zSize );
799       column.front() = fBottom->GetNode( x, y );
800       column.back()  = fTop   ->GetNode( x, y );
801     }
802   }
803
804   // compute normalized parameters of nodes on sides (PAL23189)
805   computeIJK( *fBottom, COO_X, COO_Y, /*z=*/0. );
806   computeIJK( *fRight,  COO_Y, COO_Z, /*x=*/1. );
807   computeIJK( *fTop,    COO_X, COO_Y, /*z=*/1. );
808   computeIJK( *fLeft,   COO_Y, COO_Z, /*x=*/0. );
809   computeIJK( *fFront,  COO_X, COO_Z, /*y=*/0. );
810   computeIJK( *fBack,   COO_X, COO_Z, /*y=*/1. );
811
812   StdMeshers_RenumberHelper renumHelper( aMesh, _blockRenumberHyp );
813
814   // projection points of the internal node on cube sub-shapes by which
815   // coordinates of the internal node are computed
816   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
817
818   // projections on vertices are constant
819   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
820   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
821   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
822   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
823   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
824   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
825   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
826   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
827
828   gp_XYZ params; // normalized parameters of an internal node within the unit box
829
830   if ( toRenumber )
831     for ( y = 0; y < ySize; ++y )
832     {
833       vector< const SMDS_MeshNode* >& column0y = columns[ colIndex( 0, y )];
834       for ( z = 0; z < zSize; ++z )
835         renumHelper.AddReplacingNode( column0y[ z ] );
836     }
837
838   for ( x = 1; x < xSize-1; ++x )
839   {
840     if ( toRenumber )
841     {
842       vector< const SMDS_MeshNode* >& columnX0 = columns[ colIndex( x, 0 )];
843       for ( z = 0; z < zSize; ++z )
844         renumHelper.AddReplacingNode( columnX0[ z ] );
845     }
846
847     const double rX = x / double(X);
848     for ( y = 1; y < ySize-1; ++y )
849     {
850       const double rY = y / double(Y);
851
852       // a column to fill in during z loop
853       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
854       // projection points on horizontal edges
855       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
856       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
857       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
858       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
859       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
860       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
861       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
862       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
863       // projection points on horizontal faces
864       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
865       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
866
867       if ( toRenumber )
868         renumHelper.AddReplacingNode( column[ 0 ] );
869
870       for ( z = 1; z < zSize-1; ++z ) // z loop
871       {
872         const double rZ = z / double(Z);
873
874         const gp_XYZ& pBo = fBottom->GetIJK( x, y );
875         const gp_XYZ& pTo = fTop   ->GetIJK( x, y );
876         const gp_XYZ& pFr = fFront ->GetIJK( x, z );
877         const gp_XYZ& pBa = fBack  ->GetIJK( x, z );
878         const gp_XYZ& pLe = fLeft  ->GetIJK( y, z );
879         const gp_XYZ& pRi = fRight ->GetIJK( y, z );
880         params.SetCoord( 1, 0.5 * ( pBo.X() * ( 1. - rZ ) + pTo.X() * rZ  +
881                                     pFr.X() * ( 1. - rY ) + pBa.X() * rY ));
882         params.SetCoord( 2, 0.5 * ( pBo.Y() * ( 1. - rZ ) + pTo.Y() * rZ  +
883                                     pLe.Y() * ( 1. - rX ) + pRi.Y() * rX ));
884         params.SetCoord( 3, 0.5 * ( pFr.Z() * ( 1. - rY ) + pBa.Z() * rY  +
885                                     pLe.Z() * ( 1. - rX ) + pRi.Z() * rX ));
886
887         // projection points on vertical edges
888         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );
889         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );
890         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );
891         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
892         // projection points on vertical faces
893         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );
894         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );
895         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );
896         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
897
898         // compute internal node coordinates
899         gp_XYZ coords;
900         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
901         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
902
903       } // z loop
904       if ( toRenumber )
905         renumHelper.AddReplacingNode( column[ Z ] );
906
907     } // y loop
908     if ( toRenumber )
909     {
910       vector< const SMDS_MeshNode* >& columnX0 = columns[ colIndex( x, Y )];
911       for ( z = 0; z < zSize; ++z )
912         renumHelper.AddReplacingNode( columnX0[ z ] );
913     }
914   } // x loop
915
916   if ( toRenumber )
917     for ( y = 0; y < ySize; ++y )
918     {
919       vector< const SMDS_MeshNode* >& columnXy = columns[ colIndex( X, y )];
920       for ( z = 0; z < zSize; ++z )
921         renumHelper.AddReplacingNode( columnXy[ z ] );
922     }
923
924   // side data no more needed, free memory
925   for ( int i = 0; i < 6; ++i )
926     SMESHUtils::FreeVector( aCubeSide[i]._columns );
927
928   // 5) Create hexahedrons
929   // ---------------------
930
931   for ( x = 0; x < xSize-1; ++x ) {
932     for ( y = 0; y < ySize-1; ++y ) {
933       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
934       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
935       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
936       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
937       for ( z = 0; z < zSize-1; ++z )
938       {
939         // bottom face normal of a hexa mush point outside the volume
940         if ( toRenumber )
941           helper.AddVolume(col00[z], col01[z], col01[z+1], col00[z+1],
942                            col10[z], col11[z], col11[z+1], col10[z+1]);
943         else
944           helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
945                            col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
946       }
947     }
948   }
949
950   if ( toRenumber )
951     renumHelper.DoReplaceNodes();
952
953
954   if ( _blockRenumberHyp )
955   {
956     return error( _blockRenumberHyp->CheckHypothesis( aMesh, aShape ));
957   }
958
959   return true;
960 }
961
962 //=============================================================================
963 /*!
964  *  Evaluate
965  */
966 //=============================================================================
967
968 bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
969                                   const TopoDS_Shape & aShape,
970                                   MapShapeNbElems& aResMap)
971 {
972   vector < SMESH_subMesh * >meshFaces;
973   TopTools_SequenceOfShape aFaces;
974   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
975     aFaces.Append(exp.Current());
976     SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
977     ASSERT(aSubMesh);
978     meshFaces.push_back(aSubMesh);
979   }
980   if (meshFaces.size() != 6) {
981     //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
982     static StdMeshers_CompositeHexa_3D compositeHexa(-10, aMesh.GetGen());
983     return compositeHexa.Evaluate(aMesh, aShape, aResMap);
984   }
985   
986   int i = 0;
987   for(; i<6; i++) {
988     //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
989     TopoDS_Shape aFace = aFaces.Value(i+1);
990     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
991     if( !algo ) {
992       std::vector<int> aResVec(SMDSEntity_Last);
993       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
994       SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
995       aResMap.insert(std::make_pair(sm,aResVec));
996       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
997       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
998       return false;
999     }
1000     string algoName = algo->GetName();
1001     bool isAllQuad = false;
1002     if (algoName == "Quadrangle_2D") {
1003       MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
1004       if( anIt == aResMap.end() ) continue;
1005       std::vector<int> aVec = (*anIt).second;
1006       int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1007       if( nbtri == 0 )
1008         isAllQuad = true;
1009     }
1010     if ( ! isAllQuad ) {
1011       return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
1012     }
1013   }
1014   
1015   // find number of 1d elems for 1 face
1016   int nb1d = 0;
1017   TopTools_MapOfShape Edges1;
1018   bool IsQuadratic = false;
1019   bool IsFirst = true;
1020   for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
1021     Edges1.Add(exp.Current());
1022     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
1023     if( sm ) {
1024       MapShapeNbElemsItr anIt = aResMap.find(sm);
1025       if( anIt == aResMap.end() ) continue;
1026       std::vector<int> aVec = (*anIt).second;
1027       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1028       if(IsFirst) {
1029         IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
1030         IsFirst = false;
1031       }
1032     }
1033   }
1034   // find face opposite to 1 face
1035   int OppNum = 0;
1036   for(i=2; i<=6; i++) {
1037     bool IsOpposite = true;
1038     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
1039       if( Edges1.Contains(exp.Current()) ) {
1040         IsOpposite = false;
1041         break;
1042       }
1043     }
1044     if(IsOpposite) {
1045       OppNum = i;
1046       break;
1047     }
1048   }
1049   // find number of 2d elems on side faces
1050   int nb2d = 0;
1051   for(i=2; i<=6; i++) {
1052     if( i == OppNum ) continue;
1053     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
1054     if( anIt == aResMap.end() ) continue;
1055     std::vector<int> aVec = (*anIt).second;
1056     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1057   }
1058   
1059   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
1060   std::vector<int> aVec = (*anIt).second;
1061   int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1062   int nb0d_face0 = aVec[SMDSEntity_Node];
1063
1064   std::vector<int> aResVec(SMDSEntity_Last);
1065   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1066   if(IsQuadratic) {
1067     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
1068     int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
1069     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
1070   }
1071   else {
1072     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
1073     aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
1074   }
1075   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1076   aResMap.insert(std::make_pair(sm,aResVec));
1077
1078   return true;
1079 }
1080
1081 //================================================================================
1082 /*!
1083  * \brief Computes hexahedral mesh from 2D mesh of block
1084  */
1085 //================================================================================
1086
1087 bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1088 {
1089   static StdMeshers_HexaFromSkin_3D * algo = 0;
1090   if ( !algo ) {
1091     SMESH_Gen* gen = aMesh.GetGen();
1092     algo = new StdMeshers_HexaFromSkin_3D( gen->GetANewId(), gen );
1093   }
1094   algo->InitComputeError();
1095   algo->Compute( aMesh, aHelper );
1096   return error( algo->GetComputeError());
1097 }
1098
1099 //================================================================================
1100 /*!
1101  * \brief Return true if the algorithm can mesh this shape
1102  *  \param [in] aShape - shape to check
1103  *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
1104  *              else, returns OK if at least one shape is OK
1105  */
1106 //================================================================================
1107
1108 bool StdMeshers_Hexa_3D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1109 {
1110   TopExp_Explorer exp0( aShape, TopAbs_SOLID );
1111   if ( !exp0.More() ) return false;
1112
1113   for ( ; exp0.More(); exp0.Next() )
1114   {
1115     int nbFoundShells = 0;
1116     TopExp_Explorer exp1( exp0.Current(), TopAbs_SHELL );
1117     for ( ; exp1.More(); exp1.Next(), ++nbFoundShells)
1118       if ( nbFoundShells == 2 ) break;
1119     if ( nbFoundShells != 1 ) {
1120       if ( toCheckAll ) return false;
1121       continue;
1122     }   
1123     exp1.Init( exp0.Current(), TopAbs_FACE );
1124     int nbEdges = SMESH_MesherHelper::Count( exp1.Current(), TopAbs_EDGE, /*ignoreSame=*/true );
1125     bool ok = ( nbEdges > 3 );
1126     if ( toCheckAll && !ok ) return false;
1127     if ( !toCheckAll && ok ) return true;
1128   }
1129   return toCheckAll;
1130 };
1131
1132 //=======================================================================
1133 //function : ComputePentahedralMesh
1134 //purpose  : 
1135 //=======================================================================
1136
1137 SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &          aMesh,
1138                                              const TopoDS_Shape &  aShape,
1139                                              SMESH_ProxyMesh*      proxyMesh)
1140 {
1141   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
1142   if ( proxyMesh )
1143   {
1144     err->myName = COMPERR_BAD_INPUT_MESH;
1145     err->myComment = "Can't build pentahedral mesh on viscous layers";
1146     return err;
1147   }
1148   bool bOK;
1149   StdMeshers_Penta_3D anAlgo;
1150   //
1151   bOK=anAlgo.Compute(aMesh, aShape);
1152   //
1153   err = anAlgo.GetComputeError();
1154   //
1155   if ( !bOK && anAlgo.ErrorStatus() == 5 )
1156   {
1157     static StdMeshers_Prism_3D * aPrism3D = 0;
1158     if ( !aPrism3D ) {
1159       SMESH_Gen* gen = aMesh.GetGen();
1160       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
1161     }
1162     SMESH_Hypothesis::Hypothesis_Status aStatus;
1163     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
1164       aPrism3D->InitComputeError();
1165       bOK = aPrism3D->Compute( aMesh, aShape );
1166       err = aPrism3D->GetComputeError();
1167     }
1168   }
1169   return err;
1170 }
1171
1172
1173 //=======================================================================
1174 //function : EvaluatePentahedralMesh
1175 //purpose  : 
1176 //=======================================================================
1177
1178 bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
1179                              const TopoDS_Shape & aShape,
1180                              MapShapeNbElems& aResMap)
1181 {
1182   StdMeshers_Penta_3D anAlgo;
1183   bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
1184
1185   //err = anAlgo.GetComputeError();
1186   //if ( !bOK && anAlgo.ErrorStatus() == 5 )
1187   if( !bOK ) {
1188     static StdMeshers_Prism_3D * aPrism3D = 0;
1189     if ( !aPrism3D ) {
1190       SMESH_Gen* gen = aMesh.GetGen();
1191       aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), gen );
1192     }
1193     SMESH_Hypothesis::Hypothesis_Status aStatus;
1194     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
1195       return aPrism3D->Evaluate(aMesh, aShape, aResMap);
1196     }
1197   }
1198
1199   return bOK;
1200 }