Salome HOME
Performance regression of SALOME_TESTS/Grids/smesh/bugs_18/V9
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 6227072d9ab3c5396b771dff47f6b624ffba0731..3d6249a79570ac2574e98f7b9b75965bef54b8cc 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
@@ -62,6 +62,7 @@
 #include <gp_Ax3.hxx>
 
 #include <limits>
+#include <numeric>
 
 using namespace std;
 
@@ -161,6 +162,12 @@ namespace {
     {
       return _src2tgtNodes;
     }
+    void SetEventListener( SMESH_subMesh* tgtSubMesh )
+    {
+      NSProjUtils::SetEventListener( tgtSubMesh,
+                                     _sourceHypo->GetSourceFace(),
+                                     _sourceHypo->GetSourceMesh() );
+    }
   };
   //=======================================================================
   /*!
@@ -339,7 +346,7 @@ namespace {
     // gravity center of a layer
     gp_XYZ O(0,0,0);
     int vertexCol = -1;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       O += gpXYZ( (*columns[ i ])[ z ]);
       if ( vertexCol < 0 &&
@@ -351,7 +358,7 @@ namespace {
     // Z axis
     gp_Vec Z(0,0,0);
     int iPrev = columns.size()-1;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
       gp_Vec v2( O, gpXYZ( (*columns[ i ]    )[ z ]));
@@ -363,11 +370,11 @@ namespace {
     {
       O = gpXYZ( (*columns[ vertexCol ])[ z ]);
     }
-    if ( xColumn < 0 || xColumn >= columns.size() )
+    if ( xColumn < 0 || xColumn >= (int) columns.size() )
     {
       // select a column for X dir
       double maxDist = 0;
-      for ( int i = 0; i < columns.size(); ++i )
+      for ( size_t i = 0; i < columns.size(); ++i )
       {
         double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
         if ( dist > maxDist )
@@ -454,9 +461,9 @@ namespace {
     std::advance( edgeIt, nbEdges-1 );
     TopoDS_Edge   prevE = *edgeIt;
     // bool isPrevStraight = SMESH_Algo::IsStraight( prevE );
-    int           iPrev = nbEdges - 1;
+    // int           iPrev = nbEdges - 1;
 
-    int iUnite = -1; // the first of united EDGEs
+    // int iUnite = -1; // the first of united EDGEs
 
     // analyse angles between EDGEs
     int nbCorners = 0;
@@ -524,7 +531,7 @@ namespace {
   void pointsToPython(const std::vector<gp_XYZ>& p)
   {
 #ifdef _DEBUG_
-    for ( int i = SMESH_Block::ID_V000; i < p.size(); ++i )
+    for ( size_t i = SMESH_Block::ID_V000; i < p.size(); ++i )
     {
       cout << "mesh.AddNode( " << p[i].X() << ", "<< p[i].Y() << ", "<< p[i].Z() << ") # " << i <<" " ;
       SMESH_Block::DumpShapeID( i, cout ) << endl;
@@ -560,7 +567,9 @@ StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen)
 //================================================================================
 
 StdMeshers_Prism_3D::~StdMeshers_Prism_3D()
-{}
+{
+  pointsToPython( std::vector<gp_XYZ>() ); // avoid warning: pointsToPython defined but not used
+}
 
 //=======================================================================
 //function : CheckHypothesis
@@ -933,12 +942,12 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
 
   list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
   std::list< int >::iterator     nbE = thePrism.myNbEdgesInWires.begin();
+  std::list< int > nbQuadsPerWire;
   int iE = 0;
-  double f,l;
   while ( edge != thePrism.myBottomEdges.end() )
   {
     ++iE;
-    if ( BRep_Tool::Curve( *edge, f,l ).IsNull() )
+    if ( SMESH_Algo::isDegenerated( *edge ))
     {
       edge = thePrism.myBottomEdges.erase( edge );
       --iE;
@@ -946,12 +955,14 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
     }
     else
     {
+      bool hasWallFace = false;
       TopTools_ListIteratorOfListOfShape faceIt( edgeToFaces.FindFromKey( *edge ));
       for ( ; faceIt.More(); faceIt.Next() )
       {
         const TopoDS_Face& face = TopoDS::Face( faceIt.Value() );
         if ( !thePrism.myBottom.IsSame( face ))
         {
+          hasWallFace = true;
           Prism_3D::TQuadList quadList( 1, quadAlgo->CheckNbEdges( *mesh, face ));
           if ( !quadList.back() )
             return toSM( error(TCom("Side face #") << shapeID( face )
@@ -970,12 +981,23 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
           break;
         }
       }
-      ++edge;
+      if ( hasWallFace )
+      {
+        ++edge;
+      }
+      else // seam edge (IPAL53561)
+      {
+        edge = thePrism.myBottomEdges.erase( edge );
+        --iE;
+        --(*nbE);
+      }
     }
     if ( iE == *nbE )
     {
       iE = 0;
       ++nbE;
+      int nbQuadPrev = std::accumulate( nbQuadsPerWire.begin(), nbQuadsPerWire.end(), 0 );
+      nbQuadsPerWire.push_back( thePrism.myWallQuads.size() - nbQuadPrev );
     }
   }
 
@@ -987,12 +1009,14 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   // that is not so evident in case of several WIREs in the bottom FACE
   thePrism.myRightQuadIndex.clear();
   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
-    thePrism.myRightQuadIndex.push_back( i+1 );
-  list< int >::iterator nbEinW = thePrism.myNbEdgesInWires.begin();
-  for ( int iLeft = 0; nbEinW != thePrism.myNbEdgesInWires.end(); ++nbEinW )
   {
-    thePrism.myRightQuadIndex[ iLeft + *nbEinW - 1 ] = iLeft; // 1st EDGE index of a current WIRE
-    iLeft += *nbEinW;
+    thePrism.myRightQuadIndex.push_back( i+1 ); // OK for all but the last EDGE of a WIRE
+  }
+  list< int >::iterator nbQinW = nbQuadsPerWire.begin();
+  for ( int iLeft = 0; nbQinW != nbQuadsPerWire.end(); ++nbQinW )
+  {
+    thePrism.myRightQuadIndex[ iLeft + *nbQinW - 1 ] = iLeft; // for the last EDGE of a WIRE
+    iLeft += *nbQinW;
   }
 
   while ( totalNbFaces - faceMap.Extent() > 2 )
@@ -1072,7 +1096,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   {
     // now only top and bottom FACEs are not in the faceMap
     faceMap.Add( thePrism.myBottom );
-    for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE );f.More(); f.Next() )
+    for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE ); f.More(); f.Next() )
       if ( !faceMap.Contains( f.Current() )) {
         thePrism.myTop = TopoDS::Face( f.Current() );
         break;
@@ -1110,7 +1134,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
       ( ! botSM->GetAlgo() ||
         ! _gen->Compute( *botSM->GetFather(), botSM->GetSubShape(), /*shapeOnly=*/true )))
     return error( COMPERR_BAD_INPUT_MESH,
-                  TCom( "No mesher defined to compute the face #")
+                  TCom( "No mesher defined to compute the base face #")
                   << shapeID( thePrism.myBottom ));
 
   // Make all side FACEs of thePrism meshed with quads
@@ -1196,7 +1220,8 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
     {
       const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
-      if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
+      if ( tBotNode.GetPositionType() != SMDS_TOP_FACE &&
+           myBlock.HasNodeColumn( tBotNode.myNode ))
         continue; // node is not inside the FACE
 
       // column nodes; middle part of the column are zero pointers
@@ -1286,6 +1311,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
   if ( !smDS ) return toSM( error(COMPERR_BAD_INPUT_MESH, "Null submesh"));
 
   // loop on bottom mesh faces
+  vector< const TNodeColumn* > columns;
   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
   while ( faceIt->more() )
   {
@@ -1295,24 +1321,26 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
     // find node columns for each node
     int nbNodes = face->NbCornerNodes();
-    vector< const TNodeColumn* > columns( nbNodes );
+    columns.resize( nbNodes );
     for ( int i = 0; i < nbNodes; ++i )
     {
       const SMDS_MeshNode* n = face->GetNode( i );
-      if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
+      columns[ i ] = NULL;
+
+      if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
+        columns[ i ] = myBlock.GetNodeColumn( n );
+
+      if ( !columns[ i ] )
+      {
         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
         if ( bot_column == myBotToColumnMap.end() )
-          return toSM( error(TCom("No nodes found above node ") << n->GetID() ));
-        columns[ i ] = & bot_column->second;
-      }
-      else {
-        columns[ i ] = myBlock.GetNodeColumn( n );
-        if ( !columns[ i ] )
           return toSM( error(TCom("No side nodes found above node ") << n->GetID() ));
+        columns[ i ] = & bot_column->second;
       }
     }
     // create prisms
-    AddPrisms( columns, myHelper );
+    if ( !AddPrisms( columns, myHelper ))
+      return toSM( error("Different 'vertical' discretization"));
 
   } // loop on bottom mesh faces
 
@@ -1322,7 +1350,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
   // update state of sub-meshes (mostly in order to erase improper errors)
   SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( thePrism.myShape3D );
-  SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
+  SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
   while ( smIt->more() )
   {
     sm = smIt->next();
@@ -1810,17 +1838,20 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh&         theMesh,
  */
 //================================================================================
 
-void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
+bool StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
                                      SMESH_MesherHelper*          helper)
 {
-  int nbNodes = columns.size();
-  int nbZ     = columns[0]->size();
-  if ( nbZ < 2 ) return;
+  size_t nbNodes = columns.size();
+  size_t nbZ     = columns[0]->size();
+  if ( nbZ < 2 ) return false;
+  for ( size_t i = 1; i < nbNodes; ++i )
+    if ( columns[i]->size() != nbZ )
+      return false;
 
   // find out orientation
   bool isForward = true;
   SMDS_VolumeTool vTool;
-  int z = 1;
+  size_t z = 1;
   switch ( nbNodes ) {
   case 3: {
     SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom
@@ -1906,7 +1937,7 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
     vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
     for ( z = 1; z < nbZ; ++z )
     {
-      for ( int i = 0; i < nbNodes; ++i ) {
+      for ( size_t i = 0; i < nbNodes; ++i ) {
         nodes[ i             ] = (*columns[ i ])[z+iBase1]; // bottom or top
         nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom
         // side
@@ -1920,6 +1951,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
     }
 
   } // switch ( nbNodes )
+
+  return true;
 }
 
 //================================================================================
@@ -1975,7 +2008,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
     n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap();
   }
 
-  if ( !n2nMapPtr || n2nMapPtr->size() < botSMDS->NbNodes() )
+  if ( !n2nMapPtr || (int) n2nMapPtr->size() < botSMDS->NbNodes() )
   {
     // associate top and bottom faces
     NSProjUtils::TShapeShapeMap shape2ShapeMap;
@@ -2053,7 +2086,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
   {
     const SMDS_MeshNode* botNode = bN_tN->first;
     const SMDS_MeshNode* topNode = bN_tN->second;
-    if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
+    if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
+         myBlock.HasNodeColumn( botNode ))
       continue; // wall columns are contained in myBlock
     // create node column
     Prism_3D::TNode bN( botNode );
@@ -2250,29 +2284,39 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf &             bottom
     Handle(Geom_Surface) surface = BRep_Tool::Surface( topFace, loc );
     bool isPlanar = GeomLib_IsPlanarSurface( surface ).IsPlanar();
 
-    bool isFixed = false;
     set<const SMDS_MeshNode*> fixedNodes;
-    for ( int iAttemp = 0; !isFixed && iAttemp < 10; ++iAttemp )
-    {
-      TIDSortedElemSet faces;
-      for ( faceIt = topSMDS->GetElements(); faceIt->more(); )
-        faces.insert( faces.end(), faceIt->next() );
+    TIDSortedElemSet faces;
+    for ( faceIt = topSMDS->GetElements(); faceIt->more(); )
+      faces.insert( faces.end(), faceIt->next() );
 
+    bool isOk = false;
+    for ( int isCentroidal = 0; isCentroidal < 2; ++isCentroidal )
+    {
       SMESH_MeshEditor::SmoothMethod algo =
-        iAttemp ? SMESH_MeshEditor::CENTROIDAL : SMESH_MeshEditor::LAPLACIAN;
+        isCentroidal ? SMESH_MeshEditor::CENTROIDAL : SMESH_MeshEditor::LAPLACIAN;
 
-      // smoothing
-      editor.Smooth( faces, fixedNodes, algo, /*nbIterations=*/ 10,
-                     /*theTgtAspectRatio=*/1.0, /*the2D=*/!isPlanar);
+      int nbAttempts = isCentroidal ? 1 : 10;
+      for ( int iAttemp = 0; iAttemp < nbAttempts; ++iAttemp )
+      {
+        TIDSortedElemSet workFaces = faces;
+
+        // smoothing
+        editor.Smooth( workFaces, fixedNodes, algo, /*nbIterations=*/ 10,
+                       /*theTgtAspectRatio=*/1.0, /*the2D=*/!isPlanar);
 
-      isFixed = !topHelper.IsDistorted2D( topSM, /*checkUV=*/true );
+        if (( isOk = !topHelper.IsDistorted2D( topSM, /*checkUV=*/true )) &&
+            ( !isCentroidal ))
+          break;
+      }
     }
-    if ( !isFixed )
+    if ( !isOk )
       return toSM( error( TCom("Projection from face #") << botSM->GetId()
                           << " to face #" << topSM->GetId()
                           << " failed: inverted elements created"));
   }
 
+  TProjction2dAlgo::instance( this )->SetEventListener( topSM );
+
   return true;
 }
 
@@ -2356,6 +2400,9 @@ double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePr
 
 bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism )
 {
+  if ( thePrism.myBottomEdges.size() != 4 )
+    return false;
+
   // analyse angles between edges
   double nbConcaveAng = 0, nbConvexAng = 0;
   TopoDS_Face reverseBottom = TopoDS::Face( thePrism.myBottom.Reversed() ); // see initPrism()
@@ -2411,6 +2458,8 @@ bool StdMeshers_Prism_3D::project2dMesh(const TopoDS_Face& theSrcFace,
   tgtSM->ComputeStateEngine       ( SMESH_subMesh::CHECK_COMPUTE_STATE );
   tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 
+  projector2D->SetEventListener( tgtSM );
+
   return ok;
 }
 
@@ -2490,26 +2539,34 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
   struct EdgeWithNeighbors
   {
     TopoDS_Edge _edge;
-    int         _iL, _iR;
-    EdgeWithNeighbors(const TopoDS_Edge& E, int iE, int nbE, int shift = 0 ):
-      _edge( E ),
+    int         _iBase;   /* index in a WIRE with non-base EDGEs excluded */
+    int         _iL, _iR; /* used to connect edges in a base FACE */
+    bool        _isBase;  /* is used in a base FACE */
+    EdgeWithNeighbors(const TopoDS_Edge& E, int iE, int nbE, int shift, bool isBase ):
+      _edge( E ), _iBase( iE + shift ),
       _iL( SMESH_MesherHelper::WrapIndex( iE-1, nbE ) + shift ),
-      _iR( SMESH_MesherHelper::WrapIndex( iE+1, nbE ) + shift )
+      _iR( SMESH_MesherHelper::WrapIndex( iE+1, nbE ) + shift ),
+      _isBase( isBase )
     {
     }
     EdgeWithNeighbors() {}
+    bool IsInternal() const { return !_edge.IsNull() && _edge.Orientation() == TopAbs_INTERNAL; }
   };
-  struct PrismSide
+  // PrismSide contains all FACEs linking a bottom EDGE with a top one. 
+  struct PrismSide 
   {
-    TopoDS_Face                 _face;
-    TopTools_IndexedMapOfShape *_faces; // pointer because its copy constructor is private
-    TopoDS_Edge                 _topEdge;
-    vector< EdgeWithNeighbors >*_edges;
-    int                         _iBotEdge;
-    vector< bool >              _isCheckedEdge;
+    TopoDS_Face                 _face;    // a currently treated upper FACE
+    TopTools_IndexedMapOfShape *_faces;   // all FACEs (pointer because of a private copy constructor)
+    TopoDS_Edge                 _topEdge; // a current top EDGE
+    vector< EdgeWithNeighbors >*_edges;   // all EDGEs of _face
+    int                         _iBotEdge;       // index of _topEdge within _edges
+    vector< bool >              _isCheckedEdge;  // mark EDGEs whose two owner FACEs found
     int                         _nbCheckedEdges; // nb of EDGEs whose location is defined
-    PrismSide                  *_leftSide;
+    PrismSide                  *_leftSide;       // neighbor sides
     PrismSide                  *_rightSide;
+    bool                        _isInternal; // whether this side raises from an INTERNAL EDGE
+    void SetExcluded() { _leftSide = _rightSide = NULL; }
+    bool IsExcluded() const { return !_leftSide; }
     const TopoDS_Edge& Edge( int i ) const
     {
       return (*_edges)[ i ]._edge;
@@ -2520,67 +2577,138 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
         if ( E.IsSame( Edge( i ))) return i;
       return -1;
     }
-    bool IsSideFace( const TopoDS_Shape& face ) const
+    bool IsSideFace( const TopoDS_Shape& face, const bool checkNeighbors ) const
     {
       if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
         return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
+
+      if ( checkNeighbors )
+        return (( _leftSide  && _leftSide->IsSideFace ( face, false )) ||
+                ( _rightSide && _rightSide->IsSideFace( face, false )));
+
       return false;
     }
   };
+  //--------------------------------------------------------------------------------
+  /*!
+   * \brief Return another faces sharing an edge
+   */
+  const TopoDS_Face & getAnotherFace( const TopoDS_Face& face,
+                                      const TopoDS_Edge& edge,
+                                      TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
+  {
+    TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edge ));
+    for ( ; faceIt.More(); faceIt.Next() )
+      if ( !face.IsSame( faceIt.Value() ))
+        return TopoDS::Face( faceIt.Value() );
+    return face;
+  }
+
   //--------------------------------------------------------------------------------
   /*!
    * \brief Return ordered edges of a face
    */
-  bool getEdges( const TopoDS_Face&            face,
-                 vector< EdgeWithNeighbors > & edges,
-                 const bool                    noHolesAllowed)
+  bool getEdges( const TopoDS_Face&                         face,
+                 vector< EdgeWithNeighbors > &              edges,
+                 TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge,
+                 const bool                                 noHolesAllowed)
   {
+    TopoDS_Face f = face;
+    if ( f.Orientation() != TopAbs_FORWARD &&
+         f.Orientation() != TopAbs_REVERSED )
+      f.Orientation( TopAbs_FORWARD );
     list< TopoDS_Edge > ee;
     list< int >         nbEdgesInWires;
-    int nbW = SMESH_Block::GetOrderedEdges( face, ee, nbEdgesInWires );
+    int nbW = SMESH_Block::GetOrderedEdges( f, ee, nbEdgesInWires );
     if ( nbW > 1 && noHolesAllowed )
       return false;
 
-    int iE, nbTot = 0;
-    list< TopoDS_Edge >::iterator e = ee.begin();
-    list< int >::iterator       nbE = nbEdgesInWires.begin();
+    int iE, nbTot = 0, nbBase, iBase;
+    list< TopoDS_Edge >::iterator   e = ee.begin();
+    list< int         >::iterator nbE = nbEdgesInWires.begin();
     for ( ; nbE != nbEdgesInWires.end(); ++nbE )
       for ( iE = 0; iE < *nbE; ++e, ++iE )
-        if ( SMESH_Algo::isDegenerated( *e ))
+        if ( SMESH_Algo::isDegenerated( *e )) // degenerated EDGE is never used
         {
           e = --ee.erase( e );
           --(*nbE);
           --iE;
         }
-        else
-        {
-          e->Orientation( TopAbs_FORWARD ); // for operator==() to work
-        }
 
+    vector<int> isBase;
     edges.clear();
     e = ee.begin();
     for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
     {
-      for ( iE = 0; iE < *nbE; ++e, ++iE )
-        edges.push_back( EdgeWithNeighbors( *e, iE, *nbE, nbTot ));
-      nbTot += *nbE;
+      nbBase = 0;
+      isBase.resize( *nbE );
+      list< TopoDS_Edge >::iterator eIt = e;
+      for ( iE = 0; iE < *nbE; ++eIt, ++iE )
+      {
+        isBase[ iE ] = ( getAnotherFace( face, *eIt, facesOfEdge ) != face );
+        nbBase += isBase[ iE ];
+      }
+      for ( iBase = 0, iE = 0; iE < *nbE; ++e, ++iE )
+      {
+        edges.push_back( EdgeWithNeighbors( *e, iBase, nbBase, nbTot, isBase[ iE ] ));
+        iBase += isBase[ iE ];
+      }
+      nbTot += nbBase;
+    }
+    if ( nbTot == 0 )
+      return false;
+
+    // IPAL53099. Set correct neighbors to INTERNAL EDGEs, which can be connected to
+    // EDGEs of the outer WIRE but this fact can't be detected by their order.
+    if ( nbW > 1 )
+    {
+      int iFirst = 0, iLast;
+      for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
+      {
+        iLast = iFirst + *nbE - 1;
+        TopoDS_Vertex vv[2] = { SMESH_MesherHelper::IthVertex( 0, edges[ iFirst ]._edge ),
+                                SMESH_MesherHelper::IthVertex( 1, edges[ iLast  ]._edge ) };
+        bool isConnectOk = ( vv[0].IsSame( vv[1] ));
+        if ( !isConnectOk )
+        {
+          edges[ iFirst ]._iL = edges[ iFirst ]._iBase; // connect to self
+          edges[ iLast  ]._iR = edges[ iLast ]._iBase;
+
+          // look for an EDGE of the outer WIRE connected to vv
+          TopoDS_Vertex v0, v1;
+          for ( iE = 0; iE < nbEdgesInWires.front(); ++iE )
+          {
+            v0 = SMESH_MesherHelper::IthVertex( 0, edges[ iE ]._edge );
+            v1 = SMESH_MesherHelper::IthVertex( 1, edges[ iE ]._edge );
+            if ( vv[0].IsSame( v0 ) || vv[0].IsSame( v1 ))
+              edges[ iFirst ]._iL = edges[ iE ]._iBase;
+            if ( vv[1].IsSame( v0 ) || vv[1].IsSame( v1 ))
+              edges[ iLast ]._iR = edges[ iE ]._iBase;
+          }
+        }
+        iFirst += *nbE;
+      }
     }
     return edges.size();
   }
+  
   //--------------------------------------------------------------------------------
   /*!
-   * \brief Return another faces sharing an edge
+   * \brief Return number of faces sharing given edges
    */
-  const TopoDS_Shape & getAnotherFace( const TopoDS_Face& face,
-                                       const TopoDS_Edge& edge,
-                                       TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
-  {
-    TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edge ));
-    for ( ; faceIt.More(); faceIt.Next() )
-      if ( !face.IsSame( faceIt.Value() ))
-        return faceIt.Value();
-    return face;
-  }
+  // int nbAdjacentFaces( const std::vector< EdgeWithNeighbors >&          edges,
+  //                      const TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge )
+  // {
+  //   TopTools_MapOfShape adjFaces;
+
+  //   for ( size_t i = 0; i < edges.size(); ++i )
+  //   {
+  //     TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edges[i]._edge ));
+  //     for ( ; faceIt.More(); faceIt.Next() )
+  //       adjFaces.Add( faceIt.Value() );
+  //   }
+  //   return adjFaces.Extent();
+  // }
 }
 
 //================================================================================
@@ -2603,10 +2731,10 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
     // check nb shells
     TopoDS_Shape shell;
     TopExp_Explorer shExp( sExp.Current(), TopAbs_SHELL );
-    if ( shExp.More() ) {
+    while ( shExp.More() ) {
       shell = shExp.Current();
       shExp.Next();
-      if ( shExp.More() )
+      if ( shExp.More() && BRep_Tool::IsClosed( shExp.Current() ))
         shell.Nullify();
     }
     if ( shell.IsNull() ) {
@@ -2615,7 +2743,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
     }
     // get all faces
     TopTools_IndexedMapOfShape allFaces;
-    TopExp::MapShapes( shell, TopAbs_FACE, allFaces );
+    TopExp::MapShapes( sExp.Current(), TopAbs_FACE, allFaces );
     if ( allFaces.Extent() < 3 ) {
       if ( toCheckAll ) return false;
       continue;
@@ -2632,7 +2760,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
       }
     }
 #ifdef _DEBUG_
-    TopTools_IndexedMapOfShape allShapes;
+    TopTools_IndexedMapOfShape allShapes; // usage: allShapes.FindIndex( s )
     TopExp::MapShapes( shape, allShapes );
 #endif
 
@@ -2646,32 +2774,45 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
 
     typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec;
     vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 );
-    TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ faceEdgesVec.size() ];
+    const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGEs
+    TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ];
     SMESHUtils::ArrayDeleter<TopTools_IndexedMapOfShape> delFacesOfSide( facesOfSide );
 
     // try to use each face as a bottom one
     bool prismDetected = false;
+    vector< PrismSide > sides;
     for ( int iF = 1; iF < allFaces.Extent() && !prismDetected; ++iF )
     {
       const TopoDS_Face& botF = TopoDS::Face( allFaces( iF ));
 
       TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ];
       if ( botEdges.empty() )
-        if ( !getEdges( botF, botEdges, /*noHoles=*/false ))
+        if ( !getEdges( botF, botEdges, facesOfEdge, /*noHoles=*/false ))
           break;
-      if ( allFaces.Extent()-1 <= (int) botEdges.size() )
+
+      int nbBase = 0;
+      for ( size_t iS = 0; iS < botEdges.size(); ++iS )
+        nbBase += botEdges[ iS ]._isBase;
+
+      if ( allFaces.Extent()-1 <= nbBase )
         continue; // all faces are adjacent to botF - no top FACE
 
       // init data of side FACEs
-      vector< PrismSide > sides( botEdges.size() );
-      for ( int iS = 0; iS < botEdges.size(); ++iS )
+      sides.clear();
+      sides.resize( nbBase );
+      size_t iS = 0;
+      for ( size_t iE = 0; iE < botEdges.size(); ++iE )
       {
-        sides[ iS ]._topEdge = botEdges[ iS ]._edge;
-        sides[ iS ]._face    = botF;
-        sides[ iS ]._leftSide  = & sides[ botEdges[ iS ]._iR ];
-        sides[ iS ]._rightSide = & sides[ botEdges[ iS ]._iL ];
-        sides[ iS ]._faces = & facesOfSide[ iS ];
+        if ( !botEdges[ iE ]._isBase )
+          continue;
+        sides[ iS ]._topEdge    = botEdges[ iE ]._edge;
+        sides[ iS ]._face       = botF;
+        sides[ iS ]._leftSide   = & sides[ botEdges[ iE ]._iR ];
+        sides[ iS ]._rightSide  = & sides[ botEdges[ iE ]._iL ];
+        sides[ iS ]._isInternal = botEdges[ iE ].IsInternal();
+        sides[ iS ]._faces      = & facesOfSide[ iS ];
         sides[ iS ]._faces->Clear();
+        ++iS;
       }
 
       bool isOK = true; // ok for a current botF
@@ -2699,8 +2840,9 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
                 if ( side._isCheckedEdge[ iE ] ) continue;
                 const TopoDS_Edge&      vertE = side.Edge( iE );
                 const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
-                bool isEdgeShared = adjSide->IsSideFace( neighborF );
-                if ( isEdgeShared )
+                bool isEdgeShared = (( adjSide->IsSideFace( neighborF, side._isInternal )) ||
+                                     ( adjSide == &side && neighborF.IsSame( side._face )) );
+                if ( isEdgeShared ) // vertE is shared with adjSide
                 {
                   isAdvanced = true;
                   side._isCheckedEdge[ iE ] = true;
@@ -2741,20 +2883,19 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
             {
               stop = true;
             }
-            else if ( side._leftSide != & side ) // not closed side face
+            else if ( side._leftSide != & side && // not closed side face
+                      side._leftSide->_faces->Contains( f ))
             {
-              if ( side._leftSide->_faces->Contains( f ))
-              {
-                stop = true; // probably f is the prism top face
-                side._leftSide->_face.Nullify();
-                side._leftSide->_topEdge.Nullify();
-              }
-              if ( side._rightSide->_faces->Contains( f ))
-              {
-                stop = true; // probably f is the prism top face
-                side._rightSide->_face.Nullify();
-                side._rightSide->_topEdge.Nullify();
-              }
+              stop = true; // probably f is the prism top face
+              side._leftSide->_face.Nullify();
+              side._leftSide->_topEdge.Nullify();
+            }
+            else if ( side._rightSide != & side &&
+                      side._rightSide->_faces->Contains( f ))
+            {
+              stop = true; // probably f is the prism top face
+              side._rightSide->_face.Nullify();
+              side._rightSide->_topEdge.Nullify();
             }
             if ( stop )
             {
@@ -2766,7 +2907,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
             int faceID  = allFaces.FindIndex( side._face );
             side._edges = & faceEdgesVec[ faceID ];
             if ( side._edges->empty() )
-              if ( !getEdges( side._face, * side._edges, /*noHoles=*/true ))
+              if ( !getEdges( side._face, * side._edges, facesOfEdge, /*noHoles=*/true ))
                 break;
             const int nbE = side._edges->size();
             if ( nbE >= 4 )
@@ -2779,6 +2920,10 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
               side._isCheckedEdge[ side._iBotEdge ] = true;
               side._nbCheckedEdges = 1; // bottom EDGE is known
             }
+            else // probably a triangular top face found
+            {
+              side._face.Nullify();
+            }
             side._topEdge.Nullify();
             isOK = ( !side._edges->empty() || side._faces->Extent() > 1 );
 
@@ -2811,7 +2956,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
           const TopoDS_Shape& topFace = sides[0]._faces->FindKey( nbFaces );
           size_t iS;
           for ( iS = 1; iS < sides.size(); ++iS )
-            if ( !sides[ iS ]._faces->Contains( topFace ))
+            if ( ! sides[ iS ]._faces->Contains( topFace ))
               break;
           prismDetected = ( iS == sides.size() );
         }
@@ -2937,7 +3082,6 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
   list< SMESH_subMesh* > meshedSubMesh;
   int nbFaces = 0;
   //
-  SMESH_subMesh* anyFaceSM = 0;
   SMESH_subMeshIteratorPtr smIt = mainSubMesh->getDependsOnIterator(false,true);
   while ( smIt->more() )
   {
@@ -2946,7 +3090,6 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
     if      ( face.ShapeType() > TopAbs_FACE ) break;
     else if ( face.ShapeType() < TopAbs_FACE ) continue;
     nbFaces++;
-    anyFaceSM = sm;
 
     // is quadrangle FACE?
     list< TopoDS_Edge > orderedEdges;
@@ -3076,7 +3219,7 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
 
   // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-)
   TopoDS_Vertex V000;
-  double minVal = DBL_MAX, minX, val;
+  double minVal = DBL_MAX, minX = 0, val;
   for ( TopExp_Explorer exp( botSM->GetSubShape(), TopAbs_VERTEX );
         exp.More(); exp.Next() )
   {
@@ -3219,6 +3362,9 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
+
+      if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() )
+        return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization");
     }
     // edge columns
     int id = MeshDS()->ShapeToIndex( *edgeIt );
@@ -3285,7 +3431,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
       if ( nbUnitePerEdge[ iE ] < 0 )
         continue;
       // look for already united faces
-      for ( int i = iE; i < iE + nbExraFaces; ++i )
+      for ( size_t i = iE; i < iE + nbExraFaces; ++i )
       {
         if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge
           nbExraFaces += nbUnitePerEdge[ i ];
@@ -3324,7 +3470,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     else if ( nbExraFaces > 1 ) // unite
     {
       double u0 = 0, sumLen = 0;
-      for ( int i = iE; i < iE + nbExraFaces; ++i )
+      for ( size_t i = iE; i < iE + nbExraFaces; ++i )
         sumLen += edgeLength[ i ];
 
       vector< TSideFace* >        components( nbExraFaces );
@@ -3568,7 +3714,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &
   double tol2;
   {
     Bnd_B3d bndBox;
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
       bndBox.Add( gpXYZ( columns[i]->front() ));
     tol2 = bndBox.SquareExtent() * 1e-5;
   }
@@ -3591,7 +3737,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &
     //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
 
     // check a transformation
-    for ( int i = 0; i < columns.size(); ++i )
+    for ( size_t i = 0; i < columns.size(); ++i )
     {
       gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
       gp_Pnt pz = gpXYZ( (*columns[i])[z] );
@@ -3795,7 +3941,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
   myComponents       ( other.myComponents.size() ),
   myHelper           ( *other.myHelper.GetMesh() )
 {
-  for (int i = 0 ; i < myComponents.size(); ++i )
+  for ( size_t i = 0 ; i < myComponents.size(); ++i )
     myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
 }
 
@@ -3807,7 +3953,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
 
 StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
 {
-  for (int i = 0 ; i < myComponents.size(); ++i )
+  for ( size_t i = 0 ; i < myComponents.size(); ++i )
     if ( myComponents[ i ] )
       delete myComponents[ i ];
 }
@@ -3907,7 +4053,7 @@ StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU)
   if ( myComponents.empty() )
     return const_cast<TSideFace*>( this );
 
-  int i;
+  size_t i;
   for ( i = 0; i < myComponents.size(); ++i )
     if ( U < myParams[ i ].second )
       break;
@@ -4346,9 +4492,9 @@ gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real
 void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const
 {
 #ifdef _DEBUG_
-  for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i )
+  for ( int i = 0; i < nbNodes && i < (int)myNodeColumn->size(); ++i )
     cout << (*myNodeColumn)[i]->GetID() << " ";
-  if ( nbNodes < myNodeColumn->size() )
+  if ( nbNodes < (int) myNodeColumn->size() )
     cout << myNodeColumn->back()->GetID();
 #endif
 }