Salome HOME
IPAL21120 SIGSEGV on Meshing attached Compound with Automatic Hexadralization
[modules/smesh.git] / src / StdMeshers / StdMeshers_Hexa_3D.cxx
index 9bff1bb909ed040b6dc7fcff961b6b3e8e0878de..3320a4ee94809eb43fc31c7cbed4f4527c8ace49 100644 (file)
@@ -1,38 +1,41 @@
-//  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org 
+//  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$
-
-using namespace std;
+//
 #include "StdMeshers_Hexa_3D.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"
 #include "SMESH_subMesh.hxx"
+#include "SMESH_Comment.hxx"
 
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
@@ -41,28 +44,21 @@ using namespace std;
 #include "SMDS_VolumeOfNodes.hxx"
 
 #include <TopExp.hxx>
+#include <TopExp_Explorer.hxx>
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
-#include <TopTools_ListOfShape.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
-#include <TColStd_MapOfInteger.hxx>
-
-#include <BRep_Tool.hxx>
-#include <Geom_Surface.hxx>
-#include <Geom_Curve.hxx>
-#include <Geom2d_Curve.hxx>
-#include <Handle_Geom2d_Curve.hxx>
-#include <Handle_Geom_Curve.hxx>
+#include <TopTools_ListOfShape.hxx>
+#include <TopoDS.hxx>
 #include <gp_Pnt2d.hxx>
 
 #include "utilities.h"
 #include "Utils_ExceptHandlers.hxx"
 
-//modified by NIZNHY-PKV Wed Nov 17 15:31:58 2004 f
-#include "StdMeshers_Penta_3D.hxx"
+typedef SMESH_Comment TComm;
+
+using namespace std;
 
-static bool ComputePentahedralMesh(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape);
-//modified by NIZNHY-PKV Wed Nov 17 15:32:00 2004 t
+static SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &, const TopoDS_Shape &);
 
 //=============================================================================
 /*!
@@ -70,8 +66,8 @@ static bool ComputePentahedralMesh(SMESH_Mesh & aMesh,        const TopoDS_Shape & aSha
  */
 //=============================================================================
 
-StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId,
-                                       SMESH_Gen * gen):SMESH_3D_Algo(hypId, studyId, gen)
+StdMeshers_Hexa_3D::StdMeshers_Hexa_3D(int hypId, int studyId, SMESH_Gen * gen)
+  :SMESH_3D_Algo(hypId, studyId, gen)
 {
   MESSAGE("StdMeshers_Hexa_3D::StdMeshers_Hexa_3D");
   _name = "Hexa_3D";
@@ -100,7 +96,7 @@ StdMeshers_Hexa_3D::~StdMeshers_Hexa_3D()
 bool StdMeshers_Hexa_3D::ClearAndReturn(FaceQuadStruct* theQuads[6], const bool res)
 {
   for (int i = 0; i < 6; i++) {
-    StdMeshers_Quadrangle_2D::QuadDelete(theQuads[i]);
+    delete theQuads[i];
     theQuads[i] = NULL;
   }
   return res;
@@ -114,18 +110,41 @@ bool StdMeshers_Hexa_3D::ClearAndReturn(FaceQuadStruct* theQuads[6], const bool
 //=============================================================================
 
 bool StdMeshers_Hexa_3D::CheckHypothesis
-                         (SMESH_Mesh& aMesh,
-                          const TopoDS_Shape& aShape,
+                         (SMESH_Mesh&                          aMesh,
+                          const TopoDS_Shape&                  aShape,
                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 {
-       //MESSAGE("StdMeshers_Hexa_3D::CheckHypothesis");
-
-       bool isOk = true;
-        aStatus = SMESH_Hypothesis::HYP_OK;
+  // check nb of faces in the shape
+/*  PAL16229
+  aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
+  int nbFaces = 0;
+  for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next())
+    if ( ++nbFaces > 6 )
+      break;
+  if ( nbFaces != 6 )
+    return false;
+*/
+  aStatus = SMESH_Hypothesis::HYP_OK;
+  return true;
+}
 
-       // nothing to check
+//=======================================================================
+//function : isCloser
+//purpose  : 
+//=======================================================================
 
-       return isOk;
+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;
 }
 
 //=======================================================================
@@ -135,26 +154,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->nbPts[0], quad->nbPts[2]);
-  int nbvertic = Min(quad->nbPts[1], quad->nbPts[3]);
-  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;
 }
@@ -173,16 +207,16 @@ static bool findIJ (const SMDS_MeshNode* node, const FaceQuadStruct * quad, int&
  */
 //=============================================================================
 
-bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
-       const TopoDS_Shape & aShape)throw(SALOME_Exception)
+bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
+                                 const TopoDS_Shape & aShape)// throw(SALOME_Exception)
 {
-  Unexpect aCatch(SalomeException);
+  // PAL14921. Enable catching std::bad_alloc and Standard_OutOfMemory outside
+  //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
-  //MESSAGE("---");
 
   vector < SMESH_subMesh * >meshFaces;
   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
@@ -191,12 +225,14 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
     meshFaces.push_back(aSubMesh);
   }
   if (meshFaces.size() != 6) {
-    SCRUTE(meshFaces.size());
-    return false;
+    //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)
-  //MESSAGE("---");
 
   // tool for working with quadratic elements
   SMESH_MesherHelper aTool (aMesh);
@@ -229,7 +265,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
   for (int i = 0; i < 6; i++)
     aQuads[i] = 0;
 
-  for (int i = 0; i < 6; i++) {
+  for (int i = 0; i < 6; i++)
+  {
     TopoDS_Shape aFace = meshFaces[i]->GetSubShape();
     SMESH_Algo *algo = _gen->GetAlgo(aMesh, aFace);
     string algoName = algo->GetName();
@@ -246,154 +283,88 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
       }
     }
     if ( ! isAllQuad ) {
-      //modified by NIZNHY-PKV Wed Nov 17 15:31:37 2004 f
-      bool bIsOk = ComputePentahedralMesh(aMesh, aShape);
-      return ClearAndReturn( aQuads, bIsOk );
+      SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape);
+      return ClearAndReturn( aQuads, error(err));
     }
     StdMeshers_Quadrangle_2D *quadAlgo =
       dynamic_cast < StdMeshers_Quadrangle_2D * >(algo);
     ASSERT(quadAlgo);
     try {
       aQuads[i] = quadAlgo->CheckAnd2Dcompute(aMesh, aFace, _quadraticMesh);
+      if(!aQuads[i]) {
+       return error( quadAlgo->GetComputeError());
+      }
     }
     catch(SALOME_Exception & S_ex) {
-      return ClearAndReturn( aQuads, false );
+      return ClearAndReturn( aQuads, error(COMPERR_SLM_EXCEPTION,TComm(S_ex.what()) <<
+                                           " Raised by StdMeshers_Quadrangle_2D "
+                                           " on face #" << meshDS->ShapeToIndex( aFace )));
     }
 
     // 0.2.1 - number of points on the opposite edges must be the same
-    if (aQuads[i]->nbPts[0] != aQuads[i]->nbPts[2] ||
-        aQuads[i]->nbPts[1] != aQuads[i]->nbPts[3]) {
+    if (aQuads[i]->side[0]->NbPoints() != aQuads[i]->side[2]->NbPoints() ||
+        aQuads[i]->side[1]->NbPoints() != aQuads[i]->side[3]->NbPoints()
+        /*aQuads[i]->side[0]->NbEdges() != 1 ||
+        aQuads[i]->side[1]->NbEdges() != 1 ||
+        aQuads[i]->side[2]->NbEdges() != 1 ||
+        aQuads[i]->side[3]->NbEdges() != 1*/) {
       MESSAGE("different number of points on the opposite edges of face " << i);
-      //                  ASSERT(0);
-      return ClearAndReturn( aQuads, false );
+      // Try to go into penta algorithm 'cause it has been improved.
+      SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape);
+      return ClearAndReturn( aQuads, error(err));
     }
   }
 
   // 1.  - identify faces and vertices of the "cube"
   // 1.1 - ancestor maps vertex->edges in the cube
-  //MESSAGE("---");
 
-  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
-  //MESSAGE("---");
 
   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
-  //MESSAGE("---");
-
-  int i = 0;
-  TopoDS_Edge E = aQuads[0]->edge[i];  //edge will be Y=0,Z=0 on unit cube
-  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) * (aQuads[0]->last[i] - aQuads[0]->first[i])) > 0);
-
-  if (isForward) {
-    aCube.V000 = VFirst;       // will be (0,0,0) on the unit cube
-    aCube.V100 = VLast;                // will be (1,0,0) on the unit cube
-  }
-  else {
-    aCube.V000 = VLast;
-    aCube.V100 = VFirst;
-  }
 
-  i = 1;
-  E = aQuads[0]->edge[i];
-  C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
-  TopExp::Vertices(E, VFirst, VLast);
-  isForward = (((l - f) * (aQuads[0]->last[i] - aQuads[0]->first[i])) > 0);
-  if (isForward)
-    aCube.V101 = VLast;                // will be (1,0,1) on the unit cube
-  else
-    aCube.V101 = VFirst;
-
-  i = 2;
-  E = aQuads[0]->edge[i];
-  C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
-  TopExp::Vertices(E, VFirst, VLast);
-  isForward = (((l - f) * (aQuads[0]->last[i] - aQuads[0]->first[i])) > 0);
-  if (isForward)
-    aCube.V001 = VLast;                // will be (0,0,1) on the unit cube
-  else
-    aCube.V001 = VFirst;
-
-  // 1.4 - find edge X=0, Z=0 (ancestor of V000 not in face Y=0)
-  //     - find edge X=1, Z=0 (ancestor of V100 not in face Y=0)
-  //     - find edge X=1, Z=1 (ancestor of V101 not in face Y=0) 
-  //     - find edge X=0, Z=1 (ancestor of V001 not in face Y=0)
-  //MESSAGE("---");
-
-  TopoDS_Edge E_0Y0 = EdgeNotInFace(aMesh, aShape, F, aCube.V000, MS);
-  ASSERT(!E_0Y0.IsNull());
-
-  TopoDS_Edge E_1Y0 = EdgeNotInFace(aMesh, aShape, F, aCube.V100, MS);
-  ASSERT(!E_1Y0.IsNull());
-
-  TopoDS_Edge E_1Y1 = EdgeNotInFace(aMesh, aShape, F, aCube.V101, MS);
-  ASSERT(!E_1Y1.IsNull());
-
-  TopoDS_Edge E_0Y1 = EdgeNotInFace(aMesh, aShape, F, aCube.V001, MS);
-  ASSERT(!E_0Y1.IsNull());
-
-  // 1.5 - identify the 4 vertices in face Y=1: V010, V110, V111, V011
-  //MESSAGE("---");
-
-  TopExp::Vertices(E_0Y0, VFirst, VLast);
-  if (VFirst.IsSame(aCube.V000))
-    aCube.V010 = VLast;
-  else
-    aCube.V010 = VFirst;
-
-  TopExp::Vertices(E_1Y0, VFirst, VLast);
-  if (VFirst.IsSame(aCube.V100))
-    aCube.V110 = VLast;
-  else
-    aCube.V110 = VFirst;
-
-  TopExp::Vertices(E_1Y1, VFirst, VLast);
-  if (VFirst.IsSame(aCube.V101))
-    aCube.V111 = VLast;
-  else
-    aCube.V111 = VFirst;
-
-  TopExp::Vertices(E_0Y1, VFirst, VLast);
-  if (VFirst.IsSame(aCube.V001))
-    aCube.V011 = VLast;
-  else
-    aCube.V011 = VFirst;
+  aCube.V000 = aQuads[0]->side[0]->FirstVertex(); // will be (0,0,0) on the unit cube
+  aCube.V100 = aQuads[0]->side[0]->LastVertex();  // will be (1,0,0) on the unit cube
+  aCube.V001 = aQuads[0]->side[2]->FirstVertex(); // will be (0,0,1) on the unit cube
+  aCube.V101 = aQuads[0]->side[2]->LastVertex();  // will be (1,0,1) on the unit cube
+
+  TopTools_IndexedMapOfShape MV0;
+  TopExp::MapShapes(aFace, TopAbs_VERTEX, MV0);
+
+  aCube.V010 = OppositeVertex( aCube.V000, MV0, aQuads);
+  aCube.V110 = OppositeVertex( aCube.V100, MV0, aQuads);
+  aCube.V011 = OppositeVertex( aCube.V001, MV0, aQuads);
+  aCube.V111 = OppositeVertex( aCube.V101, MV0, aQuads);
 
   // 1.6 - find remaining faces given 4 vertices
-  //MESSAGE("---");
 
   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);
-  aCube.quad_X1 = aQuads[_indX1];
 
-  //MESSAGE("---");
+  // 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
 
@@ -419,15 +390,14 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
   // 1.8 - create a 3D structure for normalized values
   
-  //MESSAGE("---");
-  int nbx = aCube.quad_Z0->nbPts[0];
-  if (cz0.a1 == 0.) nbx = aCube.quad_Z0->nbPts[1];
+  int nbx = aCube.quad_Z0->side[0]->NbPoints();
+  if (cz0.a1 == 0.) nbx = aCube.quad_Z0->side[1]->NbPoints();
  
-  int nby = aCube.quad_X0->nbPts[0];
-  if (cx0.a1 == 0.) nby = aCube.quad_X0->nbPts[1];
+  int nby = aCube.quad_X0->side[0]->NbPoints();
+  if (cx0.a1 == 0.) nby = aCube.quad_X0->side[1]->NbPoints();
  
-  int nbz = aCube.quad_Y0->nbPts[0];
-  if (cy0.a1 != 0.) nbz = aCube.quad_Y0->nbPts[1];
+  int nbz = aCube.quad_Y0->side[0]->NbPoints();
+  if (cy0.a1 != 0.) nbz = aCube.quad_Y0->side[1]->NbPoints();
 
   int i1, j1, nbxyz = nbx * nby * nbz;
   Point3DStruct *np = new Point3DStruct[nbxyz];
@@ -439,8 +409,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_X0;
     int i = 0;                         // j = x/face , k = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
 
     SMDS_NodeIteratorPtr itf= aMesh.GetSubMesh(F)->GetSubMeshDS()->GetNodes();
                        
@@ -448,7 +418,8 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -472,14 +443,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_X1;
     int i = nbx - 1;           // j = x/face , k = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
 
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -503,14 +475,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_Y0;
     int j = 0;                         // i = x/face , k = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
 
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -534,14 +507,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_Y1;
     int j = nby - 1;           // i = x/face , k = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
 
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -565,14 +539,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_Z0;
     int k = 0;                         // i = x/face , j = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
 
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -596,14 +571,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
     faceQuadStruct *quad = aCube.quad_Z1;
     int k = nbz - 1;           // i = x/face , j = y/face
-    int nbdown = quad->nbPts[0];
-    int nbright = quad->nbPts[1];
+    int nbdown = quad->side[0]->NbPoints();
+    int nbright = quad->side[1]->NbPoints();
     
     while(itf->more()) {
       const SMDS_MeshNode * node = itf->next();
       if(aTool.IsMedium(node))
         continue;
-      findIJ( node, quad, i1, j1 );
+      if ( !findIJ( node, quad, i1, j1 ))
+        return ClearAndReturn( aQuads, false );
       int ij1 = j1 * nbdown + i1;
       quad->uv_grid[ij1].node = node;
     }
@@ -740,20 +716,12 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
 
         SMDS_MeshVolume * elt;
         if ( isForw ) {
-          //elt = meshDS->AddVolume(np[n1].node, np[n2].node,
-          //                        np[n3].node, np[n4].node,
-          //                        np[n5].node, np[n6].node,
-          //                        np[n7].node, np[n8].node);
           elt = aTool.AddVolume(np[n1].node, np[n2].node,
                                 np[n3].node, np[n4].node,
                                 np[n5].node, np[n6].node,
                                 np[n7].node, np[n8].node);
         }
         else {
-          //elt = meshDS->AddVolume(np[n1].node, np[n4].node,
-          //                        np[n3].node, np[n2].node,
-          //                        np[n5].node, np[n8].node,
-          //                        np[n7].node, np[n6].node);
           elt = aTool.AddVolume(np[n1].node, np[n4].node,
                                 np[n3].node, np[n2].node,
                                 np[n5].node, np[n8].node,
@@ -765,7 +733,6 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh & aMesh,
     }
   }
   if ( np ) delete [] np;
-  //MESSAGE("End of StdMeshers_Hexa_3D::Compute()");
   return ClearAndReturn( aQuads, true );
 }
 
@@ -775,8 +742,8 @@ 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)
+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;
@@ -819,7 +786,7 @@ int StdMeshers_Hexa_3D::GetFaceIndex(SMESH_Mesh & aMesh,
                        break;
                }
        }
-       ASSERT(faceIndex > 0);
+       //IPAL21120 ASSERT(faceIndex > 0);
        //SCRUTE(faceIndex);
        return faceIndex;
 }
@@ -886,24 +853,27 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad,
        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);
-       TopoDS_Vertex VA, VB;
-       if (isForward)
-       {
-               VA = VFirst;
-               VB = VLast;
-       }
-       else
-       {
-               VA = VLast;
-               VB = VFirst;
-       }
+//     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;
+//     }
+  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))
@@ -994,8 +964,8 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad,
        conv.b2 = b2;
        conv.c2 = c2;
 
-       int nbdown = quad.nbPts[0];
-       int nbright = quad.nbPts[1];
+       int nbdown = quad.side[0]->NbPoints();
+       int nbright = quad.side[1]->NbPoints();
        conv.ia = int (a1);
        conv.ib = int (b1);
        conv.ic =
@@ -1008,48 +978,40 @@ void StdMeshers_Hexa_3D::GetConv2DCoefs(const faceQuadStruct & quad,
 //     MESSAGE("J " << conv.ja << " " << conv.jb << " " << conv.jc);
 }
 
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-ostream & StdMeshers_Hexa_3D::SaveTo(ostream & save)
-{
-  return save;
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-istream & StdMeshers_Hexa_3D::LoadFrom(istream & load)
-{
-  return load;
-}
-
-//=============================================================================
-/*!
- *  
- */
-//=============================================================================
-
-ostream & operator <<(ostream & save, StdMeshers_Hexa_3D & hyp)
-{
-  return hyp.SaveTo( save );
-}
-
-//=============================================================================
+//================================================================================
 /*!
- *  
+ * \brief Find a vertex opposite to the given vertex of aQuads[0]
+  * \param aVertex - the vertex
+  * \param aFace - the face aVertex belongs to
+  * \param aQuads - quads
+  * \retval TopoDS_Vertex - found vertex
  */
-//=============================================================================
+//================================================================================
 
-istream & operator >>(istream & load, StdMeshers_Hexa_3D & hyp)
+TopoDS_Vertex StdMeshers_Hexa_3D::OppositeVertex(const TopoDS_Vertex& aVertex,
+                                                 const TopTools_IndexedMapOfShape& aQuads0Vertices,
+                                                 FaceQuadStruct* aQuads[6])
 {
-  return hyp.LoadFrom( load );
+  int i, j;
+  for ( i = 1; i < 6; ++i )
+  {
+    TopoDS_Vertex VV[] = { aQuads[i]->side[0]->FirstVertex(),
+                           aQuads[i]->side[0]->LastVertex() , 
+                           aQuads[i]->side[2]->LastVertex() ,
+                           aQuads[i]->side[2]->FirstVertex() };
+    for ( j = 0; j < 4; ++j )
+      if ( aVertex.IsSame( VV[ j ]))
+        break;
+    if ( j < 4 ) {
+      int jPrev = j ? j - 1 : 3;
+      int jNext = (j + 1) % 4;
+      if ( aQuads0Vertices.Contains( VV[ jPrev ] ))
+        return VV[ jNext ];
+      else
+        return VV[ jPrev ];
+    }
+  }
+  return TopoDS_Vertex();
 }
 
 //modified by NIZNHY-PKV Wed Nov 17 15:34:13 2004 f
@@ -1061,26 +1023,35 @@ istream & operator >>(istream & load, StdMeshers_Hexa_3D & hyp)
 //function : ComputePentahedralMesh
 //purpose  : 
 //=======================================================================
-bool ComputePentahedralMesh(SMESH_Mesh & aMesh,        const TopoDS_Shape & aShape)
+
+SMESH_ComputeErrorPtr ComputePentahedralMesh(SMESH_Mesh &         aMesh,
+                                             const TopoDS_Shape & aShape)
 {
   //printf(" ComputePentahedralMesh HERE\n");
   //
   bool bOK;
+  SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
   //int iErr;
   StdMeshers_Penta_3D anAlgo;
   //
   bOK=anAlgo.Compute(aMesh, aShape);
-  /*
-  iErr=anAlgo.ErrorStatus();
-  
-  if (iErr) {
-    printf("  *** Error# %d\n", iErr);
-  }
-  else {
-    printf("  *** No errors# %d\n", iErr);
+  //
+  err = anAlgo.GetComputeError();
+  //
+  if ( !bOK && anAlgo.ErrorStatus() == 5 )
+  {
+    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 ) ) {
+      bOK = aPrism3D->Compute( aMesh, aShape );
+      err = aPrism3D->GetComputeError();
+    }
   }
-  */
-  return bOK;
+  return err;
 }