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