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