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