Salome HOME
Merge from V5_1_3_BR branch (07/12/09)
[modules/smesh.git] / src / StdMeshers / StdMeshers_Hexa_3D.cxx
index 33efea5e7243d429415f555ac0e3e01b818521f5..42b5488ecda855a7827d137409381fd7629e361e 100644 (file)
@@ -1,37 +1,36 @@
-//  SMESH SMESH : implementaion of SMESH idl descriptions
+//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+//  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS 
-// 
-//  This library is free software; you can redistribute it and/or 
-//  modify it under the terms of the GNU Lesser General Public 
-//  License as published by the Free Software Foundation; either 
-//  version 2.1 of the License. 
-// 
-//  This library is distributed in the hope that it will be useful, 
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of 
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
-//  Lesser General Public License for more details. 
-// 
-//  You should have received a copy of the GNU Lesser General Public 
-//  License along with this library; if not, write to the Free Software 
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA 
-// 
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License.
 //
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
 //
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
+//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+//  SMESH SMESH : implementaion of SMESH idl descriptions
 //  File   : StdMeshers_Hexa_3D.cxx
 //           Moved here from SMESH_Hexa_3D.cxx
 //  Author : Paul RASCLE, EDF
 //  Module : SMESH
-//  $Header$
-
+//
 #include "StdMeshers_Hexa_3D.hxx"
-#include "StdMeshers_Quadrangle_2D.hxx"
+#include "StdMeshers_CompositeHexa_3D.hxx"
 #include "StdMeshers_FaceSide.hxx"
 #include "StdMeshers_Penta_3D.hxx"
 #include "StdMeshers_Prism_3D.hxx"
+#include "StdMeshers_Quadrangle_2D.hxx"
 
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
@@ -49,6 +48,8 @@
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_ListOfShape.hxx>
+#include <TopTools_SequenceOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
 #include <gp_Pnt2d.hxx>
 
@@ -59,7 +60,11 @@ typedef SMESH_Comment TComm;
 
 using namespace std;
 
-static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &);
+static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &,
+                                                    const TopoDS_Shape &);
+
+static bool EvaluatePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &,
+                                    MapShapeNbElems &);
 
 //=============================================================================
 /*!
@@ -72,7 +77,7 @@ StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
 {
   MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
   _name = "Hexa_3D";
-  _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);      // 1 bit /shape type
+  _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
 }
 
 //=============================================================================
@@ -129,6 +134,25 @@ bool StdMeshers_Hexa_3D::CheckHypothesis
   return true;
 }
 
+//=======================================================================
+//function : isCloser
+//purpose  : 
+//=======================================================================
+
+inline bool isCloser(const int i, const int j, const int nbhoriz,
+                     const FaceQuadStruct* quad, const gp_Pnt2d uv,
+                     double & minDist)
+{
+  int ij = j * nbhoriz + i;
+  gp_Pnt2d uv2( quad->uv_grid[ij].u, quad->uv_grid[ij].v );
+  double dist = uv.SquareDistance( uv2 );
+  if ( dist < minDist ) {
+    minDist = dist;
+    return true;
+  }
+  return false;
+}
+
 //=======================================================================
 //function : findIJ
 //purpose  : return i,j of the node
@@ -136,26 +160,41 @@ bool StdMeshers_Hexa_3D::CheckHypothesis
 
 static bool findIJ (const SMDS_MeshNode* node, const FaceQuadStruct * quad, int& I, int& J)
 {
-  I = J = 0;
   const SMDS_FacePosition* fpos =
     static_cast<const SMDS_FacePosition*>(node->GetPosition().get());
   if ( ! fpos ) return false;
   gp_Pnt2d uv( fpos->GetUParameter(), fpos->GetVParameter() );
 
   double minDist = DBL_MAX;
-  int nbhoriz  = Min(quad->side[0]->NbPoints(), quad->side[2]->NbPoints());
-  int nbvertic = Min(quad->side[1]->NbPoints(), quad->side[3]->NbPoints());
-  for (int i = 1; i < nbhoriz - 1; i++) {
-    for (int j = 1; j < nbvertic - 1; j++) {
-      int ij = j * nbhoriz + i;
-      gp_Pnt2d uv2( quad->uv_grid[ij].u, quad->uv_grid[ij].v );
-      double dist = uv.SquareDistance( uv2 );
-      if ( dist < minDist ) {
-        minDist = dist;
-        I = i;
-        J = j;
-      }
-    }
+  const int nbhoriz  = quad->side[0]->NbPoints();
+  const int nbvertic = quad->side[1]->NbPoints();
+  I = nbhoriz/2; J = nbvertic/2;
+  int oldI, oldJ;
+  do {
+    oldI = I; oldJ = J;
+    while ( I + 2 < nbhoriz &&  isCloser( I + 1, J, nbhoriz, quad, uv, minDist ))
+      I += 1;
+    if ( I == oldI )
+      while ( I - 1 > 0     &&  isCloser( I - 1, J, nbhoriz, quad, uv, minDist ))
+        I -= 1;
+    if ( minDist < DBL_MIN )
+      break;
+
+    while ( J + 2 < nbvertic && isCloser( I, J + 1, nbhoriz, quad, uv, minDist ))
+      J += 1;
+    if ( J == oldJ )
+      while ( J - 1 > 0      && isCloser( I, J - 1, nbhoriz, quad, uv, minDist ))
+        J -= 1;
+    if ( minDist < DBL_MIN )
+      break;
+
+  } while ( I != oldI || J != oldJ );
+
+  if ( minDist > DBL_MIN ) {
+    for (int i = 1; i < nbhoriz - 1; i++)
+      for (int j = 1; j < nbvertic - 1; j++)
+        if ( isCloser( i, j, nbhoriz, quad, uv, minDist ))
+          I = i, J = j;
   }
   return true;
 }
@@ -181,7 +220,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   //Unexpect aCatch(SalomeException);
   MESSAGE("StdMeshers_Hexa_3D::Compute");
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
-  
+
   // 0.  - shape and face mesh verification
   // 0.1 - shape must be a solid (or a shell) with 6 faces
 
@@ -191,8 +230,13 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     ASSERT(aSubMesh);
     meshFaces.push_back(aSubMesh);
   }
-  if (meshFaces.size() != 6)
-    return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
+  if (meshFaces.size() != 6) {
+    //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
+    static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
+    if ( !compositeHexa.Compute( aMesh, aShape ))
+      return error( compositeHexa.GetComputeError() );
+    return true;
+  }
 
   // 0.2 - is each face meshed with Quadrangle_2D? (so, with a wire of 4 edges)
 
@@ -254,7 +298,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     try {
       aQuads[i] = quadAlgo->CheckAnd2Dcompute(aMesh, aFace, _quadraticMesh);
       if(!aQuads[i]) {
-       return error( quadAlgo->GetComputeError());
+        return error( quadAlgo->GetComputeError());
       }
     }
     catch(SALOME_Exception & S_ex) {
@@ -280,13 +324,13 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   // 1.  - identify faces and vertices of the "cube"
   // 1.1 - ancestor maps vertex->edges in the cube
 
-  TopTools_IndexedDataMapOfShapeListOfShape MS;
-  TopExp::MapShapesAndAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, MS);
+//   TopTools_IndexedDataMapOfShapeListOfShape MS;
+//   TopExp::MapShapesAndAncestors(aShape, TopAbs_VERTEX, TopAbs_EDGE, MS);
 
   // 1.2 - first face is choosen as face Y=0 of the unit cube
 
   const TopoDS_Shape & aFace = meshFaces[0]->GetSubShape();
-  const TopoDS_Face & F = TopoDS::Face(aFace);
+  //const TopoDS_Face & F = TopoDS::Face(aFace);
 
   // 1.3 - identify the 4 vertices of the face Y=0: V000, V100, V101, V001
 
@@ -296,7 +340,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   aCube.V101 = aQuads[0]->side[2]->LastVertex();  // will be (1,0,1) on the unit cube
 
   TopTools_IndexedMapOfShape MV0;
-  TopExp::MapShapes(F, TopAbs_VERTEX, MV0);
+  TopExp::MapShapes(aFace, TopAbs_VERTEX, MV0);
 
   aCube.V010 = OppositeVertex( aCube.V000, MV0, aQuads);
   aCube.V110 = OppositeVertex( aCube.V100, MV0, aQuads);
@@ -306,32 +350,32 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   // 1.6 - find remaining faces given 4 vertices
 
   int _indY0 = 0;
-  aCube.quad_Y0 = aQuads[_indY0];
-
   int _indY1 = GetFaceIndex(aMesh, aShape, meshFaces,
                             aCube.V010, aCube.V011, aCube.V110, aCube.V111);
-  aCube.quad_Y1 = aQuads[_indY1];
-
   int _indZ0 = GetFaceIndex(aMesh, aShape, meshFaces,
                             aCube.V000, aCube.V010, aCube.V100, aCube.V110);
-  aCube.quad_Z0 = aQuads[_indZ0];
-
   int _indZ1 = GetFaceIndex(aMesh, aShape, meshFaces,
                             aCube.V001, aCube.V011, aCube.V101, aCube.V111);
-  aCube.quad_Z1 = aQuads[_indZ1];
-
   int _indX0 = GetFaceIndex(aMesh, aShape, meshFaces,
                             aCube.V000, aCube.V001, aCube.V010, aCube.V011);
-  aCube.quad_X0 = aQuads[_indX0];
-
   int _indX1 = GetFaceIndex(aMesh, aShape, meshFaces,
                             aCube.V100, aCube.V101, aCube.V110, aCube.V111);
+
+  // IPAL21120: SIGSEGV on Meshing attached Compound with Automatic Hexadralization
+  if ( _indY1 < 1 || _indZ0 < 1 || _indZ1 < 1 || _indX0 < 1 || _indX1 < 1 )
+    return error(COMPERR_BAD_SHAPE);
+
+  aCube.quad_Y0 = aQuads[_indY0];
+  aCube.quad_Y1 = aQuads[_indY1];
+  aCube.quad_Z0 = aQuads[_indZ0];
+  aCube.quad_Z1 = aQuads[_indZ1];
+  aCube.quad_X0 = aQuads[_indX0];
   aCube.quad_X1 = aQuads[_indX1];
 
   // 1.7 - get convertion coefs from face 2D normalized to 3D normalized
 
-  Conv2DStruct cx0;                    // for face X=0
-  Conv2DStruct cx1;                    // for face X=1
+  Conv2DStruct cx0;                     // for face X=0
+  Conv2DStruct cx1;                     // for face X=1
   Conv2DStruct cy0;
   Conv2DStruct cy1;
   Conv2DStruct cz0;
@@ -370,12 +414,12 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     const TopoDS_Face & F = TopoDS::Face(meshFaces[_indX0]->GetSubShape());
 
     faceQuadStruct *quad = aCube.quad_X0;
-    int i = 0;                         // j = x/face , k = y/face
+    int i = 0;                          // j = x/face , k = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
 
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
-                       
+                        
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
@@ -389,8 +433,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int j = cx0.ia * i1 + cx0.ib * j1 + cx0.ic;    // j = x/face
-        int k = cx0.ja * i1 + cx0.jb * j1 + cx0.jc;    // k = y/face
+        int j = cx0.ia * i1 + cx0.ib * j1 + cx0.ic;     // j = x/face
+        int k = cx0.ja * i1 + cx0.jb * j1 + cx0.jc;     // k = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -404,7 +448,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
 
     faceQuadStruct *quad = aCube.quad_X1;
-    int i = nbx - 1;           // j = x/face , k = y/face
+    int i = nbx - 1;            // j = x/face , k = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
 
@@ -421,8 +465,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int j = cx1.ia * i1 + cx1.ib * j1 + cx1.ic;    // j = x/face
-        int k = cx1.ja * i1 + cx1.jb * j1 + cx1.jc;    // k = y/face
+        int j = cx1.ia * i1 + cx1.ib * j1 + cx1.ic;     // j = x/face
+        int k = cx1.ja * i1 + cx1.jb * j1 + cx1.jc;     // k = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -436,7 +480,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
 
     faceQuadStruct *quad = aCube.quad_Y0;
-    int j = 0;                         // i = x/face , k = y/face
+    int j = 0;                          // i = x/face , k = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
 
@@ -453,8 +497,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int i = cy0.ia * i1 + cy0.ib * j1 + cy0.ic;    // i = x/face
-        int k = cy0.ja * i1 + cy0.jb * j1 + cy0.jc;    // k = y/face
+        int i = cy0.ia * i1 + cy0.ib * j1 + cy0.ic;     // i = x/face
+        int k = cy0.ja * i1 + cy0.jb * j1 + cy0.jc;     // k = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -468,7 +512,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
 
     faceQuadStruct *quad = aCube.quad_Y1;
-    int j = nby - 1;           // i = x/face , k = y/face
+    int j = nby - 1;            // i = x/face , k = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
 
@@ -485,8 +529,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int i = cy1.ia * i1 + cy1.ib * j1 + cy1.ic;    // i = x/face
-        int k = cy1.ja * i1 + cy1.jb * j1 + cy1.jc;    // k = y/face
+        int i = cy1.ia * i1 + cy1.ib * j1 + cy1.ic;     // i = x/face
+        int k = cy1.ja * i1 + cy1.jb * j1 + cy1.jc;     // k = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -500,7 +544,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
 
     faceQuadStruct *quad = aCube.quad_Z0;
-    int k = 0;                         // i = x/face , j = y/face
+    int k = 0;                          // i = x/face , j = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
 
@@ -517,8 +561,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int i = cz0.ia * i1 + cz0.ib * j1 + cz0.ic;    // i = x/face
-        int j = cz0.ja * i1 + cz0.jb * j1 + cz0.jc;    // j = y/face
+        int i = cz0.ia * i1 + cz0.ib * j1 + cz0.ic;     // i = x/face
+        int j = cz0.ja * i1 + cz0.jb * j1 + cz0.jc;     // j = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -532,7 +576,7 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
 
     faceQuadStruct *quad = aCube.quad_Z1;
-    int k = nbz - 1;           // i = x/face , j = y/face
+    int k = nbz - 1;            // i = x/face , j = y/face
     int nbdown = quad->side[0]->NbPoints();
     int nbright = quad->side[1]->NbPoints();
     
@@ -549,8 +593,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     for (int i1 = 0; i1 < nbdown; i1++)
       for (int j1 = 0; j1 < nbright; j1++) {
         int ij1 = j1 * nbdown + i1;
-        int i = cz1.ia * i1 + cz1.ib * j1 + cz1.ic;    // i = x/face
-        int j = cz1.ja * i1 + cz1.jb * j1 + cz1.jc;    // j = y/face
+        int i = cz1.ia * i1 + cz1.ib * j1 + cz1.ic;     // i = x/face
+        int j = cz1.ja * i1 + cz1.jb * j1 + cz1.jc;     // j = y/face
         int ijk = k * nbx * nby + j * nbx + i;
         //MESSAGE(" "<<ij1<<" "<<i<<" "<<j<<" "<<ijk);
         np[ijk].node = quad->uv_grid[ij1].node;
@@ -611,9 +655,9 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
         GetPoint(p1yz, nbx - 1, j, k, nbx, nby, nbz, np, meshDS);
 
         int ijk = k * nbx * nby + j * nbx + i;
-        double x = double (i) / double (nbx - 1);      // *** seulement
-        double y = double (j) / double (nby - 1);      // *** maillage
-        double z = double (k) / double (nbz - 1);      // *** regulier
+        double x = double (i) / double (nbx - 1);       // *** seulement
+        double y = double (j) / double (nby - 1);       // *** maillage
+        double z = double (k) / double (nbz - 1);       // *** regulier
 
         Pt3 X;
         for (int i = 0; i < 3; i++) {
@@ -698,6 +742,127 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   return ClearAndReturn( aQuads, true );
 }
 
+
+//=============================================================================
+/*!
+ *  Evaluate
+ */
+//=============================================================================
+
+bool StdMeshers_Hexa_3D::Evaluate(SMESH_Mesh & aMesh,
+                                  const TopoDS_Shape & aShape,
+                                  MapShapeNbElems& aResMap)
+{
+  vector < SMESH_subMesh * >meshFaces;
+  TopTools_SequenceOfShape aFaces;
+  for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
+    aFaces.Append(exp.Current());
+    SMESH_subMesh *aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
+    ASSERT(aSubMesh);
+    meshFaces.push_back(aSubMesh);
+  }
+  if (meshFaces.size() != 6) {
+    //return error(COMPERR_BAD_SHAPE, TComm(meshFaces.size())<<" instead of 6 faces in a block");
+    static StdMeshers_CompositeHexa_3D compositeHexa(-10, 0, aMesh.GetGen());
+    return compositeHexa.Evaluate(aMesh, aShape, aResMap);
+  }
+  
+  int i = 0;
+  for(; i<6; i++) {
+    //TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
+    TopoDS_Shape aFace = aFaces.Value(i+1);
+    SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
+    if( !algo ) {
+      std::vector<int> aResVec(SMDSEntity_Last);
+      for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+      SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+      aResMap.insert(std::make_pair(sm,aResVec));
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+      return false;
+    }
+    string algoName = algo->GetName();
+    bool isAllQuad = false;
+    if (algoName == "Quadrangle_2D") {
+      MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
+      if( anIt == aResMap.end() ) continue;
+      std::vector<int> aVec = (*anIt).second;
+      int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
+      if( nbtri == 0 )
+        isAllQuad = true;
+    }
+    if ( ! isAllQuad ) {
+      return EvaluatePentahedralMesh(aMesh, aShape, aResMap);
+    }
+  }
+  
+  // find number of 1d elems for 1 face
+  int nb1d = 0;
+  TopTools_MapOfShape Edges1;
+  bool IsQuadratic = false;
+  bool IsFirst = true;
+  for (TopExp_Explorer exp(aFaces.Value(1), TopAbs_EDGE); exp.More(); exp.Next()) {
+    Edges1.Add(exp.Current());
+    SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
+    if( sm ) {
+      MapShapeNbElemsItr anIt = aResMap.find(sm);
+      if( anIt == aResMap.end() ) continue;
+      std::vector<int> aVec = (*anIt).second;
+      nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
+      if(IsFirst) {
+        IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
+        IsFirst = false;
+      }
+    }
+  }
+  // find face opposite to 1 face
+  int OppNum = 0;
+  for(i=2; i<=6; i++) {
+    bool IsOpposite = true;
+    for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
+      if( Edges1.Contains(exp.Current()) ) {
+        IsOpposite = false;
+        break;
+      }
+    }
+    if(IsOpposite) {
+      OppNum = i;
+      break;
+    }
+  }
+  // find number of 2d elems on side faces
+  int nb2d = 0;
+  for(i=2; i<=6; i++) {
+    if( i == OppNum ) continue;
+    MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
+    if( anIt == aResMap.end() ) continue;
+    std::vector<int> aVec = (*anIt).second;
+    nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+  }
+  
+  MapShapeNbElemsItr anIt = aResMap.find( meshFaces[0] );
+  std::vector<int> aVec = (*anIt).second;
+  int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
+  int nb0d_face0 = aVec[SMDSEntity_Node];
+
+  std::vector<int> aResVec(SMDSEntity_Last);
+  for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
+  if(IsQuadratic) {
+    aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0 * ( nb2d/nb1d );
+    int nb1d_face0_int = ( nb2d_face0*4 - nb1d ) / 2;
+    aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
+  }
+  else {
+    aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
+    aResVec[SMDSEntity_Hexa] = nb2d_face0 * ( nb2d/nb1d );
+  }
+  SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
+  aResMap.insert(std::make_pair(sm,aResVec));
+
+  return true;
+}
+
+
 //=============================================================================
 /*!
  *  
@@ -707,12 +872,12 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
 void StdMeshers_Hexa_3D::GetPoint(Pt3 p, int i, int j, int k, int nbx, int nby, int nbz,
                                   Point3DStruct * np, const SMESHDS_Mesh * meshDS)
 {
-       int ijk = k * nbx * nby + j * nbx + i;
-       const SMDS_MeshNode * node = np[ijk].node;
-       p[0] = node->X();
-       p[1] = node->Y();
-       p[2] = node->Z();
-       //MESSAGE(" "<<i<<" "<<j<<" "<<k<<" "<<p[0]<<" "<<p[1]<<" "<<p[2]);
+        int ijk = k * nbx * nby + j * nbx + i;
+        const SMDS_MeshNode * node = np[ijk].node;
+        p[0] = node->X();
+        p[1] = node->Y();
+        p[2] = node->Z();
+        //MESSAGE(" "<<i<<" "<<j<<" "<<k<<" "<<p[0]<<" "<<p[1]<<" "<<p[2]);
 }
 
 //=============================================================================
@@ -722,35 +887,35 @@ void StdMeshers_Hexa_3D::GetPoint(Pt3 p, int i, int j, int k, int nbx, int nby,
 //=============================================================================
 
 int StdMeshers_Hexa_3D::GetFaceIndex(SMESH_Mesh & aMesh,
-       const TopoDS_Shape & aShape,
-       const vector < SMESH_subMesh * >&meshFaces,
-       const TopoDS_Vertex & V0,
-       const TopoDS_Vertex & V1,
-       const TopoDS_Vertex & V2, const TopoDS_Vertex & V3)
+        const TopoDS_Shape & aShape,
+        const vector < SMESH_subMesh * >&meshFaces,
+        const TopoDS_Vertex & V0,
+        const TopoDS_Vertex & V1,
+        const TopoDS_Vertex & V2, const TopoDS_Vertex & V3)
 {
-       //MESSAGE("StdMeshers_Hexa_3D::GetFaceIndex");
-       int faceIndex = -1;
-       for (int i = 1; i < 6; i++)
-       {
-               const TopoDS_Shape & aFace = meshFaces[i]->GetSubShape();
-               //const TopoDS_Face& F = TopoDS::Face(aFace);
-               TopTools_IndexedMapOfShape M;
-               TopExp::MapShapes(aFace, TopAbs_VERTEX, M);
-               bool verticesInShape = false;
-               if (M.Contains(V0))
-                       if (M.Contains(V1))
-                               if (M.Contains(V2))
-                                       if (M.Contains(V3))
-                                               verticesInShape = true;
-               if (verticesInShape)
-               {
-                       faceIndex = i;
-                       break;
-               }
-       }
-       ASSERT(faceIndex > 0);
-       //SCRUTE(faceIndex);
-       return faceIndex;
+        //MESSAGE("StdMeshers_Hexa_3D::GetFaceIndex");
+        int faceIndex = -1;
+        for (int i = 1; i < 6; i++)
+        {
+                const TopoDS_Shape & aFace = meshFaces[i]->GetSubShape();
+                //const TopoDS_Face& F = TopoDS::Face(aFace);
+                TopTools_IndexedMapOfShape M;
+                TopExp::MapShapes(aFace, TopAbs_VERTEX, M);
+                bool verticesInShape = false;
+                if (M.Contains(V0))
+                        if (M.Contains(V1))
+                                if (M.Contains(V2))
+                                        if (M.Contains(V3))
+                                                verticesInShape = true;
+                if (verticesInShape)
+                {
+                        faceIndex = i;
+                        break;
+                }
+        }
+        //IPAL21120 ASSERT(faceIndex > 0);
+        //SCRUTE(faceIndex);
+        return faceIndex;
 }
 
 //=============================================================================
@@ -760,46 +925,46 @@ int StdMeshers_Hexa_3D::GetFaceIndex(SMESH_Mesh & aMesh,
 //=============================================================================
 
 TopoDS_Edge
-       StdMeshers_Hexa_3D::EdgeNotInFace(SMESH_Mesh & aMesh,
-       const TopoDS_Shape & aShape,
-       const TopoDS_Face & aFace,
-       const TopoDS_Vertex & aVertex,
-       const TopTools_IndexedDataMapOfShapeListOfShape & MS)
+        StdMeshers_Hexa_3D::EdgeNotInFace(SMESH_Mesh & aMesh,
+        const TopoDS_Shape & aShape,
+        const TopoDS_Face & aFace,
+        const TopoDS_Vertex & aVertex,
+        const TopTools_IndexedDataMapOfShapeListOfShape & MS)
 {
-       //MESSAGE("StdMeshers_Hexa_3D::EdgeNotInFace");
-       TopTools_IndexedDataMapOfShapeListOfShape MF;
-       TopExp::MapShapesAndAncestors(aFace, TopAbs_VERTEX, TopAbs_EDGE, MF);
-       const TopTools_ListOfShape & ancestorsInSolid = MS.FindFromKey(aVertex);
-       const TopTools_ListOfShape & ancestorsInFace = MF.FindFromKey(aVertex);
-//     SCRUTE(ancestorsInSolid.Extent());
-//     SCRUTE(ancestorsInFace.Extent());
-       ASSERT(ancestorsInSolid.Extent() == 6); // 6 (edges doublees)
-       ASSERT(ancestorsInFace.Extent() == 2);
-
-       TopoDS_Edge E;
-       E.Nullify();
-       TopTools_ListIteratorOfListOfShape its(ancestorsInSolid);
-       for (; its.More(); its.Next())
-       {
-               TopoDS_Shape ancestor = its.Value();
-               TopTools_ListIteratorOfListOfShape itf(ancestorsInFace);
-               bool isInFace = false;
-               for (; itf.More(); itf.Next())
-               {
-                       TopoDS_Shape ancestorInFace = itf.Value();
-                       if (ancestorInFace.IsSame(ancestor))
-                       {
-                               isInFace = true;
-                               break;
-                       }
-               }
-               if (!isInFace)
-               {
-                       E = TopoDS::Edge(ancestor);
-                       break;
-               }
-       }
-       return E;
+        //MESSAGE("StdMeshers_Hexa_3D::EdgeNotInFace");
+        TopTools_IndexedDataMapOfShapeListOfShape MF;
+        TopExp::MapShapesAndAncestors(aFace, TopAbs_VERTEX, TopAbs_EDGE, MF);
+        const TopTools_ListOfShape & ancestorsInSolid = MS.FindFromKey(aVertex);
+        const TopTools_ListOfShape & ancestorsInFace = MF.FindFromKey(aVertex);
+//      SCRUTE(ancestorsInSolid.Extent());
+//      SCRUTE(ancestorsInFace.Extent());
+        ASSERT(ancestorsInSolid.Extent() == 6); // 6 (edges doublees)
+        ASSERT(ancestorsInFace.Extent() == 2);
+
+        TopoDS_Edge E;
+        E.Nullify();
+        TopTools_ListIteratorOfListOfShape its(ancestorsInSolid);
+        for (; its.More(); its.Next())
+        {
+                TopoDS_Shape ancestor = its.Value();
+                TopTools_ListIteratorOfListOfShape itf(ancestorsInFace);
+                bool isInFace = false;
+                for (; itf.More(); itf.Next())
+                {
+                        TopoDS_Shape ancestorInFace = itf.Value();
+                        if (ancestorInFace.IsSame(ancestor))
+                        {
+                                isInFace = true;
+                                break;
+                        }
+                }
+                if (!isInFace)
+                {
+                        E = TopoDS::Edge(ancestor);
+                        break;
+                }
+        }
+        return E;
 }
 
 //=============================================================================
@@ -809,135 +974,135 @@ TopoDS_Edge
 //=============================================================================
 
 void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad,
-       const TopoDS_Shape & aShape,
-       const TopoDS_Vertex & V0,
-       const TopoDS_Vertex & V1,
-       const TopoDS_Vertex & V2, const TopoDS_Vertex & V3, Conv2DStruct & conv)
+        const TopoDS_Shape & aShape,
+        const TopoDS_Vertex & V0,
+        const TopoDS_Vertex & V1,
+        const TopoDS_Vertex & V2, const TopoDS_Vertex & V3, Conv2DStruct & conv)
 {
-//     MESSAGE("StdMeshers_Hexa_3D::GetConv2DCoefs");
-//     const TopoDS_Face & F = TopoDS::Face(aShape);
-//     TopoDS_Edge E = quad.edge[0];
-//     double f, l;
-//     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
-//     TopoDS_Vertex VFirst, VLast;
-//     TopExp::Vertices(E, VFirst, VLast);     // corresponds to f and l
-//     bool isForward = (((l - f) * (quad.last[0] - quad.first[0])) > 0);
+//      MESSAGE("StdMeshers_Hexa_3D::GetConv2DCoefs");
+//      const TopoDS_Face & F = TopoDS::Face(aShape);
+//      TopoDS_Edge E = quad.edge[0];
+//      double f, l;
+//      Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
+//      TopoDS_Vertex VFirst, VLast;
+//      TopExp::Vertices(E, VFirst, VLast);     // corresponds to f and l
+//      bool isForward = (((l - f) * (quad.last[0] - quad.first[0])) > 0);
   TopoDS_Vertex VA, VB;
-//     if (isForward)
-//     {
-//             VA = VFirst;
-//             VB = VLast;
-//     }
-//     else
-//     {
-//             VA = VLast;
-//             VB = VFirst;
-//     }
+//      if (isForward)
+//      {
+//              VA = VFirst;
+//              VB = VLast;
+//      }
+//      else
+//      {
+//              VA = VLast;
+//              VB = VFirst;
+//      }
   VA = quad.side[0]->FirstVertex();
   VB = quad.side[0]->LastVertex();
   
-       int a1, b1, c1, a2, b2, c2;
-       if (VA.IsSame(V0))
-               if (VB.IsSame(V1))
-               {
-                       a1 = 1;
-                       b1 = 0;
-                       c1 = 0;                         // x
-                       a2 = 0;
-                       b2 = 1;
-                       c2 = 0;                         // y
-               }
-               else
-               {
-                       ASSERT(VB.IsSame(V3));
-                       a1 = 0;
-                       b1 = 1;
-                       c1 = 0;                         // y
-                       a2 = 1;
-                       b2 = 0;
-                       c2 = 0;                         // x
-               }
-       if (VA.IsSame(V1))
-               if (VB.IsSame(V2))
-               {
-                       a1 = 0;
-                       b1 = -1;
-                       c1 = 1;                         // 1-y
-                       a2 = 1;
-                       b2 = 0;
-                       c2 = 0;                         // x
-               }
-               else
-               {
-                       ASSERT(VB.IsSame(V0));
-                       a1 = -1;
-                       b1 = 0;
-                       c1 = 1;                         // 1-x
-                       a2 = 0;
-                       b2 = 1;
-                       c2 = 0;                         // y
-               }
-       if (VA.IsSame(V2))
-               if (VB.IsSame(V3))
-               {
-                       a1 = -1;
-                       b1 = 0;
-                       c1 = 1;                         // 1-x
-                       a2 = 0;
-                       b2 = -1;
-                       c2 = 1;                         // 1-y
-               }
-               else
-               {
-                       ASSERT(VB.IsSame(V1));
-                       a1 = 0;
-                       b1 = -1;
-                       c1 = 1;                         // 1-y
-                       a2 = -1;
-                       b2 = 0;
-                       c2 = 1;                         // 1-x
-               }
-       if (VA.IsSame(V3))
-               if (VB.IsSame(V0))
-               {
-                       a1 = 0;
-                       b1 = 1;
-                       c1 = 0;                         // y
-                       a2 = -1;
-                       b2 = 0;
-                       c2 = 1;                         // 1-x
-               }
-               else
-               {
-                       ASSERT(VB.IsSame(V2));
-                       a1 = 1;
-                       b1 = 0;
-                       c1 = 0;                         // x
-                       a2 = 0;
-                       b2 = -1;
-                       c2 = 1;                         // 1-y
-               }
-//     MESSAGE("X = " << c1 << "+ " << a1 << "*x + " << b1 << "*y");
-//     MESSAGE("Y = " << c2 << "+ " << a2 << "*x + " << b2 << "*y");
-       conv.a1 = a1;
-       conv.b1 = b1;
-       conv.c1 = c1;
-       conv.a2 = a2;
-       conv.b2 = b2;
-       conv.c2 = c2;
-
-       int nbdown = quad.side[0]->NbPoints();
-       int nbright = quad.side[1]->NbPoints();
-       conv.ia = int (a1);
-       conv.ib = int (b1);
-       conv.ic =
-               int (c1 * a1 * a1) * (nbdown - 1) + int (c1 * b1 * b1) * (nbright - 1);
-       conv.ja = int (a2);
-       conv.jb = int (b2);
-       conv.jc =
-               int (c2 * a2 * a2) * (nbdown - 1) + int (c2 * b2 * b2) * (nbright - 1);
-//     MESSAGE("I " << conv.ia << " " << conv.ib << " " << conv.ic);
-//     MESSAGE("J " << conv.ja << " " << conv.jb << " " << conv.jc);
+        int a1, b1, c1, a2, b2, c2;
+        if (VA.IsSame(V0))
+                if (VB.IsSame(V1))
+                {
+                        a1 = 1;
+                        b1 = 0;
+                        c1 = 0;                         // x
+                        a2 = 0;
+                        b2 = 1;
+                        c2 = 0;                         // y
+                }
+                else
+                {
+                        ASSERT(VB.IsSame(V3));
+                        a1 = 0;
+                        b1 = 1;
+                        c1 = 0;                         // y
+                        a2 = 1;
+                        b2 = 0;
+                        c2 = 0;                         // x
+                }
+        if (VA.IsSame(V1))
+                if (VB.IsSame(V2))
+                {
+                        a1 = 0;
+                        b1 = -1;
+                        c1 = 1;                         // 1-y
+                        a2 = 1;
+                        b2 = 0;
+                        c2 = 0;                         // x
+                }
+                else
+                {
+                        ASSERT(VB.IsSame(V0));
+                        a1 = -1;
+                        b1 = 0;
+                        c1 = 1;                         // 1-x
+                        a2 = 0;
+                        b2 = 1;
+                        c2 = 0;                         // y
+                }
+        if (VA.IsSame(V2))
+                if (VB.IsSame(V3))
+                {
+                        a1 = -1;
+                        b1 = 0;
+                        c1 = 1;                         // 1-x
+                        a2 = 0;
+                        b2 = -1;
+                        c2 = 1;                         // 1-y
+                }
+                else
+                {
+                        ASSERT(VB.IsSame(V1));
+                        a1 = 0;
+                        b1 = -1;
+                        c1 = 1;                         // 1-y
+                        a2 = -1;
+                        b2 = 0;
+                        c2 = 1;                         // 1-x
+                }
+        if (VA.IsSame(V3))
+                if (VB.IsSame(V0))
+                {
+                        a1 = 0;
+                        b1 = 1;
+                        c1 = 0;                         // y
+                        a2 = -1;
+                        b2 = 0;
+                        c2 = 1;                         // 1-x
+                }
+                else
+                {
+                        ASSERT(VB.IsSame(V2));
+                        a1 = 1;
+                        b1 = 0;
+                        c1 = 0;                         // x
+                        a2 = 0;
+                        b2 = -1;
+                        c2 = 1;                         // 1-y
+                }
+//      MESSAGE("X = " << c1 << "+ " << a1 << "*x + " << b1 << "*y");
+//      MESSAGE("Y = " << c2 << "+ " << a2 << "*x + " << b2 << "*y");
+        conv.a1 = a1;
+        conv.b1 = b1;
+        conv.c1 = c1;
+        conv.a2 = a2;
+        conv.b2 = b2;
+        conv.c2 = c2;
+
+        int nbdown = quad.side[0]->NbPoints();
+        int nbright = quad.side[1]->NbPoints();
+        conv.ia = int (a1);
+        conv.ib = int (b1);
+        conv.ic =
+                int (c1 * a1 * a1) * (nbdown - 1) + int (c1 * b1 * b1) * (nbright - 1);
+        conv.ja = int (a2);
+        conv.jb = int (b2);
+        conv.jc =
+                int (c2 * a2 * a2) * (nbdown - 1) + int (c2 * b2 * b2) * (nbright - 1);
+//      MESSAGE("I " << conv.ia << " " << conv.ib << " " << conv.ic);
+//      MESSAGE("J " << conv.ja << " " << conv.jb << " " << conv.jc);
 }
 
 //================================================================================
@@ -1009,6 +1174,7 @@ SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &         aMesh,
     }
     SMESH_Hypothesis::Hypothesis_Status aStatus;
     if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
+      aPrism3D->InitComputeError();
       bOK = aPrism3D->Compute( aMesh, aShape );
       err = aPrism3D->GetComputeError();
     }
@@ -1017,3 +1183,33 @@ SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &         aMesh,
 }
 
 
+//=======================================================================
+//function : EvaluatePentahedralMesh
+//purpose  : 
+//=======================================================================
+
+bool EvaluatePentahedralMesh(SMESH_Mesh & aMesh,
+                             const TopoDS_Shape & aShape,
+                             MapShapeNbElems& aResMap)
+{
+  StdMeshers_Penta_3D anAlgo;
+  bool bOK = anAlgo.Evaluate(aMesh, aShape, aResMap);
+
+  //err = anAlgo.GetComputeError();
+  //if ( !bOK && anAlgo.ErrorStatus() == 5 )
+  if( !bOK ) {
+    static StdMeshers_Prism_3D * aPrism3D = 0;
+    if ( !aPrism3D ) {
+      SMESH_Gen* gen = aMesh.GetGen();
+      aPrism3D = new StdMeshers_Prism_3D( gen->GetANewId(), 0, gen );
+    }
+    SMESH_Hypothesis::Hypothesis_Status aStatus;
+    if ( aPrism3D->CheckHypothesis( aMesh, aShape, aStatus ) ) {
+      return aPrism3D->Evaluate(aMesh, aShape, aResMap);
+    }
+  }
+
+  return bOK;
+}
+
+