Salome HOME
IPAL52499: Prismatic mesh is not computed on a prismatic shape
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 8b7347d7b9702cc27cd8b02b523e7da47039e776..5e253829ab1af52a83a2b8a25489a086df85af4a 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
@@ -51,6 +51,7 @@
 #include <Geom2d_Line.hxx>
 #include <GeomLib_IsPlanarSurface.hxx>
 #include <Geom_Curve.hxx>
+#include <TColStd_DataMapOfIntegerInteger.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
@@ -162,6 +163,12 @@ namespace {
     {
       return _src2tgtNodes;
     }
+    void SetEventListener( SMESH_subMesh* tgtSubMesh )
+    {
+      NSProjUtils::SetEventListener( tgtSubMesh,
+                                     _sourceHypo->GetSourceFace(),
+                                     _sourceHypo->GetSourceMesh() );
+    }
   };
   //=======================================================================
   /*!
@@ -561,7 +568,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
@@ -864,7 +873,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
 
 
     // TODO. there are other ways to find out the source FACE:
-    // propagation, topological similarity, ect.
+    // propagation, topological similarity, etc...
 
     // simply try to mesh all not meshed SOLIDs
     if ( meshedFaces.empty() )
@@ -936,11 +945,10 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   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;
@@ -948,12 +956,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 )
@@ -972,7 +982,16 @@ 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 )
     {
@@ -1116,7 +1135,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
@@ -1155,7 +1174,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
   // Projections on the top and bottom faces are taken from nodes existing
   // on these faces; find correspondence between bottom and top nodes
-  myUseBlock = false;
+  myUseBlock = false; // is set to true if projection is done using "block approach"
   myBotToColumnMap.clear();
   if ( !assocOrProjBottom2Top( bottomToTopTrsf, thePrism ) ) // it also fills myBotToColumnMap
     return false;
@@ -1163,38 +1182,60 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
   // Create nodes inside the block
 
-  // use transformation (issue 0020680, IPAL0052499)
-  StdMeshers_Sweeper sweeper;
-  double tol;
-  bool allowHighBndError;
-
   if ( !myUseBlock )
   {
+    // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach
+    StdMeshers_Sweeper sweeper;
+
     // load boundary nodes into sweeper
     bool dummy;
+    const SMDS_MeshNode* prevN0 = 0, *prevN1 = 0;
     list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin();
     for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
     {
       int edgeID = meshDS->ShapeToIndex( *edge );
       TParam2ColumnMap* u2col = const_cast<TParam2ColumnMap*>
         ( myBlock.GetParam2ColumnMap( edgeID, dummy ));
-      TParam2ColumnMap::iterator u2colIt = u2col->begin();
-      for ( ; u2colIt != u2col->end(); ++u2colIt )
+
+      TParam2ColumnMap::iterator u2colIt = u2col->begin(), u2colEnd = u2col->end();
+      const SMDS_MeshNode* n0 = u2colIt->second[0];
+      const SMDS_MeshNode* n1 = u2col->rbegin()->second[0];
+      if ( n0 == prevN0 || n0 == prevN1 ) ++u2colIt;
+      if ( n1 == prevN0 || n1 == prevN1 ) --u2colEnd;
+      prevN0 = n0; prevN1 = n1;
+
+      for ( ; u2colIt != u2colEnd; ++u2colIt )
         sweeper.myBndColumns.push_back( & u2colIt->second );
     }
-    // load node columns inside the bottom face
+    // load node columns inside the bottom FACE
     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
+    sweeper.myIntColumns.reserve( myBotToColumnMap.size() );
     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
       sweeper.myIntColumns.push_back( & bot_column->second );
 
-    tol = getSweepTolerance( thePrism );
-    allowHighBndError = !isSimpleBottom( thePrism );
-  }
+    myHelper->SetElementsOnShape( true );
 
-  if ( !myUseBlock && sweeper.ComputeNodes( *myHelper, tol, allowHighBndError ))
-  {
+    // If all "vertical" EDGEs are straight, then all nodes of an internal node column
+    // are located on a line connecting the top node and the bottom node.
+    bool isStrightColunm = allVerticalEdgesStraight( thePrism );
+    if ( !isStrightColunm )
+    {
+      double tol = getSweepTolerance( thePrism );
+      bool allowHighBndError = !isSimpleBottom( thePrism );
+      myUseBlock = !sweeper.ComputeNodes( *myHelper, tol, allowHighBndError );
+    }
+    else if ( sweeper.CheckSameZ() )
+    {
+      myUseBlock = !sweeper.ComputeNodesOnStraightSameZ( *myHelper );
+    }
+    else
+    {
+      myUseBlock = !sweeper.ComputeNodesOnStraight( *myHelper, thePrism.myBottom, thePrism.myTop );
+    }
+    myHelper->SetElementsOnShape( false );
   }
-  else // use block approach
+
+  if ( myUseBlock ) // use block approach
   {
     // loop on nodes inside the bottom face
     Prism_3D::TNode prevBNode;
@@ -1202,7 +1243,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
@@ -1306,16 +1348,17 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
     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
@@ -1330,7 +1373,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();
@@ -1583,7 +1626,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
       if ( ! fSM->IsMeshComputed() )
       {
         // Top EDGEs must be projections from the bottom ones
-        // to compute stuctured quad mesh on wall FACEs
+        // to compute structured quad mesh on wall FACEs
         // ---------------------------------------------------
         const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge(0);
         const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ].grid->Edge(0);
@@ -2066,7 +2109,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 );
@@ -2263,29 +2307,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;
+
+      int nbAttempts = isCentroidal ? 1 : 10;
+      for ( int iAttemp = 0; iAttemp < nbAttempts; ++iAttemp )
+      {
+        TIDSortedElemSet workFaces = faces;
 
-      // smoothing
-      editor.Smooth( faces, fixedNodes, algo, /*nbIterations=*/ 10,
-                     /*theTgtAspectRatio=*/1.0, /*the2D=*/!isPlanar);
+        // 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;
 }
 
@@ -2369,6 +2423,9 @@ double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePr
 
 bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism )
 {
+  if ( thePrism.myNbEdgesInWires.front() != 4 )
+    return false;
+
   // analyse angles between edges
   double nbConcaveAng = 0, nbConvexAng = 0;
   TopoDS_Face reverseBottom = TopoDS::Face( thePrism.myBottom.Reversed() ); // see initPrism()
@@ -2398,6 +2455,43 @@ bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism )
   return true;
 }
 
+//=======================================================================
+//function : allVerticalEdgesStraight
+//purpose  : Defines if all "vertical" EDGEs are straight
+//=======================================================================
+
+bool StdMeshers_Prism_3D::allVerticalEdgesStraight( const Prism_3D::TPrismTopo& thePrism )
+{
+  for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
+  {
+    const Prism_3D::TQuadList& quads = thePrism.myWallQuads[i];
+    Prism_3D::TQuadList::const_iterator quadIt = quads.begin();
+    TopoDS_Edge prevQuadEdge;
+    for ( ; quadIt != quads.end(); ++quadIt )
+    {
+      StdMeshers_FaceSidePtr rightSide = (*quadIt)->side[ QUAD_RIGHT_SIDE ];
+
+      if ( !prevQuadEdge.IsNull() &&
+           !SMESH_Algo::IsContinuous( rightSide->Edge( 0 ), prevQuadEdge ))
+        return false;
+
+      for ( int iE = 0; iE < rightSide->NbEdges(); ++iE )
+      {
+        const TopoDS_Edge & rightE = rightSide->Edge( iE );
+        if ( !SMESH_Algo::IsStraight( rightE, /*degenResult=*/true ))
+          return false;
+
+        if ( iE > 0 &&
+             !SMESH_Algo::IsContinuous( rightSide->Edge( iE-1 ), rightE ))
+          return false;
+
+        prevQuadEdge = rightE;
+      }
+    }
+  }
+  return true;
+}
+
 //=======================================================================
 //function : project2dMesh
 //purpose  : Project mesh faces from a source FACE of one prism (theSrcFace)
@@ -2424,6 +2518,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;
 }
 
@@ -2503,14 +2599,18 @@ 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 ),
-      _iL( SMESH_MesherHelper::WrapIndex( iE-1, nbE ) + shift ),
-      _iR( SMESH_MesherHelper::WrapIndex( iE+1, nbE ) + shift )
+    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, Max( 1, nbE )) + shift ),
+      _iR( SMESH_MesherHelper::WrapIndex( iE+1, Max( 1, nbE )) + shift ),
+      _isBase( isBase )
     {
     }
     EdgeWithNeighbors() {}
+    bool IsInternal() const { return !_edge.IsNull() && _edge.Orientation() == TopAbs_INTERNAL; }
   };
   // PrismSide contains all FACEs linking a bottom EDGE with a top one. 
   struct PrismSide 
@@ -2524,6 +2624,7 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
     int                         _nbCheckedEdges; // nb of EDGEs whose location is defined
     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
@@ -2536,85 +2637,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 WIREs connected to vv
+          TopoDS_Vertex v0, v1;
+          for ( iE = 0; iE < iFirst; ++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
-   */
-  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 number of faces sharing given edges
    */
-  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();
-  }
+  // 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();
+  // }
 }
 
 //================================================================================
@@ -2637,10 +2791,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() ) {
@@ -2649,7 +2803,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;
@@ -2666,7 +2820,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
 
@@ -2693,22 +2847,32 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
 
       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
       sides.clear();
-      sides.resize( botEdges.size() );
-      for ( size_t iS = 0; iS < botEdges.size(); ++iS )
+      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
@@ -2736,8 +2900,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 )               // vertE is shared with adjSide
+                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;
@@ -2778,20 +2943,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 )
             {
@@ -2803,7 +2967,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 )
@@ -2978,7 +3142,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() )
   {
@@ -2987,7 +3150,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;
@@ -3117,7 +3279,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() )
   {
@@ -3156,9 +3318,11 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
                     "Non-quadrilateral faces are not opposite"));
 
     // check that the found top and bottom FACEs are opposite
+    TopTools_IndexedMapOfShape topEdgesMap( thePrism.myBottomEdges.size() );
+    TopExp::MapShapes( thePrism.myTop, topEdgesMap );
     list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
     for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
-      if ( myHelper->IsSubShape( *edge, thePrism.myTop ))
+      if ( topEdgesMap.Contains( *edge ))
         return toSM( error
                      (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
                       "Non-quadrilateral faces are not opposite"));
@@ -3260,6 +3424,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 );
@@ -4625,7 +4792,7 @@ void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints,
 
 //================================================================================
 /*!
- * \brief Creates internal nodes of the prism
+ * \brief Create internal nodes of the prism
  */
 //================================================================================
 
@@ -4787,8 +4954,8 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
     }
   }
 
-  //centerIntErrorIsSmall = true;
-  //bndErrorIsSmall = true;
+  centerIntErrorIsSmall = true; // 3D_mesh_Extrusion_00/A3
+  bndErrorIsSmall = true;
   if ( !centerIntErrorIsSmall )
   {
     // Compensate the central error; continue adding projection
@@ -4905,3 +5072,257 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
 
   return true;
 }
+
+//================================================================================
+/*!
+ * \brief Check if all nodes of each layers have same logical Z
+ */
+//================================================================================
+
+bool StdMeshers_Sweeper::CheckSameZ()
+{
+  myZColumns.resize( myBndColumns.size() );
+  fillZColumn( myZColumns[0], *myBndColumns[0] );
+
+  bool sameZ = true;
+  const double tol = 0.1 * 1./ myBndColumns[0]->size();
+
+  // check columns based on VERTEXes
+
+  vector< int > vertexIndex;
+  vertexIndex.push_back( 0 );
+  for ( size_t iC = 1; iC < myBndColumns.size() &&  sameZ; ++iC )
+  {
+    if ( myBndColumns[iC]->front()->GetPosition()->GetDim() > 0 )
+      continue; // not on VERTEX
+
+    vertexIndex.push_back( iC );
+    fillZColumn( myZColumns[iC], *myBndColumns[iC] );
+
+    for ( size_t iZ = 0; iZ < myZColumns[0].size() &&  sameZ; ++iZ )
+      sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol );
+  }
+
+  // check columns based on EDGEs, one per EDGE
+
+  for ( size_t i = 1; i < vertexIndex.size() &&  sameZ; ++i )
+  {
+    if ( vertexIndex[i] - vertexIndex[i-1] < 2 )
+      continue;
+
+    int iC = ( vertexIndex[i] + vertexIndex[i-1] ) / 2;
+    fillZColumn( myZColumns[iC], *myBndColumns[iC] );
+
+    for ( size_t iZ = 0; iZ < myZColumns[0].size() &&  sameZ; ++iZ )
+      sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol );
+  }
+
+  if ( sameZ )
+  {
+    myZColumns.resize(1);
+  }
+  else
+  {
+    for ( size_t iC = 1; iC < myBndColumns.size(); ++iC )
+      fillZColumn( myZColumns[iC], *myBndColumns[iC] );
+  }
+
+  return sameZ;
+}
+
+//================================================================================
+/*!
+ * \brief Create internal nodes of the prism all located on straight lines with
+ *        the same distribution along the lines.
+ */
+//================================================================================
+
+bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper )
+{
+  TZColumn& z = myZColumns[0];
+
+  for ( size_t i = 0; i < myIntColumns.size(); ++i )
+  {
+    TNodeColumn& nodes = *myIntColumns[i];
+    SMESH_NodeXYZ n0( nodes[0] ), n1( nodes.back() );
+
+    for ( size_t iZ = 0; iZ < z.size(); ++iZ )
+    {
+      gp_XYZ p = n0 * ( 1 - z[iZ] ) + n1 * z[iZ];
+      nodes[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() );
+    }
+  }
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Create internal nodes of the prism all located on straight lines with
+ *        different distributions along the lines.
+ */
+//================================================================================
+
+bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper,
+                                                 const TopoDS_Face&  botFace,
+                                                 const TopoDS_Face&  topFace )
+{
+  // get data to create a Morph
+  UVPtStructVec botUV( myBndColumns.size() + 1 );
+  UVPtStructVec topUV( myBndColumns.size() + 1 );
+  for ( size_t i = 0; i < myBndColumns.size(); ++i )
+  {
+    TNodeColumn& nodes = *myBndColumns[i];
+    botUV[i].node = nodes[0];
+    botUV[i].SetUV( helper.GetNodeUV( botFace, nodes[0] ));
+    topUV[i].node = nodes.back();
+    topUV[i].SetUV( helper.GetNodeUV( topFace, nodes.back() ));
+    botUV[i].node->setIsMarked( true );
+  }
+  botUV.back() = botUV[0];
+  topUV.back() = topUV[0];
+
+  TopoDS_Edge dummyE;
+  TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, botFace, dummyE, helper.GetMesh() ));
+  TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, topFace, dummyE, helper.GetMesh() ));
+
+  // use Morph to make delauney mesh on the FACEs. Locating of a node within a
+  // delauney triangle will be used to get a weighted Z.
+  NSProjUtils::Morph botDelauney( botWires );
+  NSProjUtils::Morph topDelauney( topWires );
+
+  if ( helper.GetIsQuadratic() )
+  {
+    // mark all medium nodes of faces on botFace to avoid their treating
+    SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( botFace );
+    SMDS_ElemIteratorPtr eIt = smDS->GetElements();
+    while ( eIt->more() )
+    {
+      const SMDS_MeshElement* e = eIt->next();
+      for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i )
+        e->GetNode( i )->setIsMarked( true );
+    }
+  }
+
+  // map to get a node column by a bottom node
+  TColStd_DataMapOfIntegerInteger iNode2iCol( myIntColumns.size() );
+
+  // un-mark nodes to treat (internal bottom nodes); later we will mark treated nodes
+  for ( size_t i = 0; i < myIntColumns.size(); ++i )
+  {
+    const SMDS_MeshNode* botNode = myIntColumns[i]->front();
+    botNode->setIsMarked( false );
+    iNode2iCol.Bind( botNode->GetID(), i );
+  }
+
+  const int botFaceID = helper.GetMesh()->GetSubMesh( botFace )->GetId();
+  const SMDS_MeshNode     *botNode, *topNode;
+  const BRepMesh_Triangle *botTria, *topTria;
+  double botBC[3], topBC[3]; // barycentric coordinates
+  int    botTriaNodes[3], topTriaNodes[3];
+  bool   checkUV = true;
+
+  // a queue of bottom nodes with starting delauney triangles
+  NSProjUtils::Morph::TNodeTriaList botNoTriQueue;
+
+  size_t iBndN = 1; // index of a bottom boundary node
+  int nbNodesToProcess = myIntColumns.size();
+  while ( nbNodesToProcess > 0 )
+  {
+    while ( !botNoTriQueue.empty() ) // treat all nodes in the queue
+    {
+      botNode = botNoTriQueue.front().first;
+      botTria = botNoTriQueue.front().second;
+      botNoTriQueue.pop_front();
+      if ( botNode->isMarked() )
+        continue;
+      --nbNodesToProcess;
+      botNode->setIsMarked( true );
+
+      TNodeColumn* column = myIntColumns[ iNode2iCol( botNode->GetID() )];
+
+      // find a delauney triangle containing the botNode
+      gp_XY botUV = helper.GetNodeUV( botFace, botNode, NULL, &checkUV );
+      botUV  *= botDelauney.GetScale();
+      botTria = botDelauney.FindTriangle( botUV, botTria, botBC, botTriaNodes );
+      if ( !botTria )
+        return false;
+
+      // find a delauney triangle containing the topNode
+      topNode = column->back();
+      gp_XY topUV = helper.GetNodeUV( topFace, topNode, NULL, &checkUV );
+      topUV *= topDelauney.GetScale();
+      // get a starting triangle basing on that top and bot boundary nodes have same index
+      topTria = topDelauney.GetTriangleNear( botTriaNodes[0] );
+      topTria = topDelauney.FindTriangle( topUV, topTria, topBC, topTriaNodes );
+      if ( !topTria )
+        return false;
+
+      // create nodes along a line
+      SMESH_NodeXYZ botP( botNode ), topP( topNode);
+      for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ )
+      {
+        // use barycentric coordinates as weight of Z of boundary columns
+        double botZ = 0, topZ = 0;
+        for ( int i = 0; i < 3; ++i )
+        {
+          botZ += botBC[i] * myZColumns[ botTriaNodes[i]-1 ][ iZ ];
+          topZ += topBC[i] * myZColumns[ topTriaNodes[i]-1 ][ iZ ];
+        }
+        double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 );
+        double z = botZ * ( 1 - rZ ) + topZ * rZ;
+        gp_XYZ p = botP * ( 1 - z  ) + topP * z;
+        (*column)[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() );
+      }
+
+      // add neighbor nodes to the queue
+      botDelauney.AddCloseNodes( botNode, botTria, botFaceID, botNoTriQueue );
+    }
+
+    if ( nbNodesToProcess > 0 ) // fill the queue
+    {
+      // assure that all bot nodes are visited
+      for ( ; iBndN-1 < myBndColumns.size() &&  botNoTriQueue.empty();  ++iBndN )
+      {
+        botTria = botDelauney.GetTriangleNear( iBndN );
+        const SMDS_MeshNode*  bndNode = botDelauney.GetBndNodes()[ iBndN ];
+        botDelauney.AddCloseNodes( bndNode, botTria, botFaceID, botNoTriQueue );
+      }
+      if ( botNoTriQueue.empty() )
+      {
+        for ( size_t i = 0; i < myIntColumns.size(); ++i )
+        {
+          botNode = myIntColumns[i]->front();
+          if ( !botNode->isMarked() )
+            botNoTriQueue.push_back( make_pair( botNode, botTria ));
+        }
+      }
+    }
+  }
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Compute Z of nodes of a straight column
+ */
+//================================================================================
+
+void StdMeshers_Sweeper::fillZColumn( TZColumn&    zColumn,
+                                      TNodeColumn& nodes )
+{
+  if ( zColumn.size() == nodes.size() - 2 )
+    return;
+
+  gp_Pnt p0 = SMESH_NodeXYZ( nodes[0] );
+  gp_Vec line( p0, SMESH_NodeXYZ( nodes.back() ));
+  double len2 = line.SquareMagnitude();
+
+  zColumn.resize( nodes.size() - 2 );
+  for ( size_t i = 0; i < zColumn.size(); ++i )
+  {
+    gp_Vec vec( p0, SMESH_NodeXYZ( nodes[ i+1] ));
+    zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line
+  }
+}