Salome HOME
merge V5_1_4
[modules/smesh.git] / src / StdMeshers / StdMeshers_HexaFromSkin_3D.cxx
index 28b628c5c3d476f98272092ea12d5be06effe1a0..6e8448e0a2092d0bf0226dc03dc77a8393f55df2 100644 (file)
@@ -51,7 +51,7 @@ namespace
 {
   enum EBoxSides //!< sides of the block
     {
-      B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED
+      B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
     };
   const char* SBoxSides[] = //!< names of block sides
     {
@@ -59,7 +59,7 @@ namespace
     };
   enum EQuadEdge //!< edges of quadrangle side
     {
-      Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED
+      Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
     };
 
 
@@ -84,9 +84,42 @@ namespace
     return true;
   }
 
+  //================================================================================
+  /*!
+   * \brief Description of node used to detect corner nodes
+   */
+  struct _NodeDescriptor
+  {
+    int nbInverseFaces, nbNodesInInverseFaces;
+
+    _NodeDescriptor(const SMDS_MeshNode* n): nbInverseFaces(0), nbNodesInInverseFaces(0)
+    {
+      if ( n )
+      {
+        set<const SMDS_MeshNode*> nodes;
+        SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
+        while ( fIt->more() )
+        {
+          const SMDS_MeshElement* face = fIt->next();
+          nodes.insert( face->begin_nodes(), face->end_nodes() );
+          nbInverseFaces++;
+        }
+        nbNodesInInverseFaces = nodes.size();
+      }
+    }
+    bool operator==(const _NodeDescriptor& other) const
+    {
+      return
+        nbInverseFaces == other.nbInverseFaces &&
+        nbNodesInInverseFaces == other.nbNodesInInverseFaces;
+    }
+  };
+
   //================================================================================
   /*!
    * \brief return true if a node is at block corner
+   *
+   * This check is valid for simple cases only
    */
   //================================================================================
 
@@ -103,7 +136,7 @@ namespace
 
   bool isQuadrangle(const SMDS_MeshElement* e)
   {
-    return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 ));
+    return ( e && e->NbCornerNodes() == 4 );
   }
 
   //================================================================================
@@ -189,15 +222,15 @@ namespace
     int                          _nbBlocksExpected;
     int                          _nbBlocksFound;
 
-#ifdef _DEBUG_
-#define _grid_access_(args) _grid.at( args )
+#ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
+#define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
 #else
-#define _grid_access_(args) _grid[ args ]
+#define _grid_access_(pobj, i) pobj->_grid[ i ]
 #endif
     //!< Return node at XY
-    const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); }
+    const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
     //!< Set node at XY
-    void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; }
+    void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
     //!< Return an edge
     SMESH_OrientedLink getEdge(EQuadEdge edge) const
     {
@@ -242,7 +275,7 @@ namespace
     //!< return coordinates by XY
     gp_XYZ xyz(int x, int y) const
     {
-      return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) );
+      return SMESH_MeshEditor::TNodeXYZ( _grid_access_(_side, _index( x, y )) );
     }
     //!< safely return a node by XY
     const SMDS_MeshNode* node(int x, int y) const
@@ -259,7 +292,7 @@ namespace
     //!< Return a corner node
     const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
     {
-      return _side->_grid_access_( _index.corner( isXMax, isYMax ));
+      return _grid_access_(_side, _index.corner( isXMax, isYMax ));
     }
     //!< return its size in nodes
     int getHoriSize() const { return _index.xSize(); }
@@ -276,14 +309,29 @@ namespace
    */
   struct _Block
   {
-    _OrientedBlockSide _side[6]; // 6 sides of a sub-block
+    _OrientedBlockSide        _side[6]; // 6 sides of a sub-block
+    set<const SMDS_MeshNode*> _corners;
 
     const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
+    bool setSide( int i, const _OrientedBlockSide& s)
+    {
+      if (( _side[i] = s ))
+      {
+        _corners.insert( s.cornerNode(0,0));
+        _corners.insert( s.cornerNode(1,0));
+        _corners.insert( s.cornerNode(0,1));
+        _corners.insert( s.cornerNode(1,1));
+      }
+      return s;
+    }
+    void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
     bool hasSide( const _OrientedBlockSide& s) const
     {
       if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
       return false;
     }
+    int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
+    bool isValid() const;
   };
   //================================================================================
   /*!
@@ -300,7 +348,9 @@ namespace
     const SMESH_Comment& error() const { return _error; }
 
   private:
-    bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad);
+    bool fillSide( _BlockSide&             side,
+                   const SMDS_MeshElement* cornerQuad,
+                   const SMDS_MeshNode*    cornerNode);
     bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
                              const SMDS_MeshNode*    n1,
                              const SMDS_MeshNode*    n2,
@@ -309,16 +359,14 @@ namespace
                              bool alongN1N2 );
     _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
                                       EQuadEdge sharedSideEdge1,
-                                      EQuadEdge sharedSideEdge2);
+                                      EQuadEdge sharedSideEdge2,
+                                      bool      withGeometricAnalysis);
     //!< update own data and data of the side bound to block
     void setSideBoundToBlock( _BlockSide& side )
     {
-      side._nbBlocksFound++;
-      if ( side.isBound() )
-      {
-        for ( int e = 0; e < int(Q_UNDEFINED); ++e )
+      if ( side._nbBlocksFound++, side.isBound() )
+        for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
           _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
-      }
     }
     //!< store reason of error
     int error(const SMESH_Comment& reason) { _error = reason; return 0; }
@@ -383,7 +431,7 @@ namespace
           _allSides.push_back( _BlockSide() );
 
         _BlockSide& side = _allSides.back();
-        if ( !fillSide( side, face ) )
+        if ( !fillSide( side, face, *corner ) )
         {
           if ( !_error.empty() )
             return false;
@@ -394,13 +442,10 @@ namespace
             for ( int isYMax = 0; isYMax < 2; ++isYMax )
             {
               const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
-              if ( !isCornerNode( nCorner ))
-                return BAD_MESH_ERR;
-              //_corner2sides[ nCorner ].insert( &side );
               corners.push_back( nCorner );
               cornerFaces.insert( side.getCornerFace( nCorner ));
             }
-          for ( int e = 0; e < int(Q_UNDEFINED); ++e )
+          for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
             _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
 
           nbFacesOnSides += side.getNbFaces();
@@ -425,23 +470,33 @@ namespace
       return BAD_MESH_ERR;
     if ( _allSides.back()._grid.empty() )
       _allSides.pop_back();
+    _DUMP_("Nb detected sides "<< _allSides.size());
 
     // ---------------------------
     // Organize sides into blocks
     // ---------------------------
 
-    // analyse sharing of sides by blocks
-    int nbBlockSides = 0; // nb of block sides taking into account their sharing
-    list < _BlockSide >::iterator sideIt = _allSides.begin();
-    for ( ; sideIt != _allSides.end(); ++sideIt )
+    // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
+    int nbBlockSides = 0; // total nb of block sides taking into account their sharing
+    multimap<int, _BlockSide* > sortedSides;
     {
-      _BlockSide& side = *sideIt;
-      bool isSharedSide = true;
-      for ( int e = 0; e < int(Q_UNDEFINED) && isSharedSide; ++e )
-        isSharedSide = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size() > 2;
-      side._nbBlocksFound = 0;
-      side._nbBlocksExpected = isSharedSide ? 2 : 1;
-      nbBlockSides += side._nbBlocksExpected;
+      list < _BlockSide >::iterator sideIt = _allSides.begin();
+      for ( ; sideIt != _allSides.end(); ++sideIt )
+      {
+        _BlockSide& side = *sideIt;
+        bool isSharedSide = true;
+        int nbAdjacent = 0;
+        for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
+        {
+          int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
+          nbAdjacent += nbAdj;
+          isSharedSide = ( nbAdj > 2 );
+        }
+        side._nbBlocksFound = 0;
+        side._nbBlocksExpected = isSharedSide ? 2 : 1;
+        nbBlockSides += side._nbBlocksExpected;
+        sortedSides.insert( make_pair( nbAdjacent, & side ));
+      }
     }
 
     // find sides of each block
@@ -449,9 +504,9 @@ namespace
     while ( nbBlockSides >= 6 )
     {
       // get any side not bound to all blocks it belongs to
-      sideIt = _allSides.begin();
-      while ( sideIt != _allSides.end() && sideIt->isBound())
-        ++sideIt;
+      multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
+      while ( i_side != sortedSides.end() && i_side->second->isBound())
+        ++i_side;
 
       // start searching for block sides from the got side
       bool ok = true;
@@ -459,37 +514,40 @@ namespace
         _blocks.resize( _blocks.size() + 1 );
 
       _Block& block = _blocks.back();
-      block._side[B_FRONT] = &(*sideIt);
-      setSideBoundToBlock( *sideIt );
+      block.setSide( B_FRONT, i_side->second );
+      setSideBoundToBlock( *i_side->second );
       nbBlockSides--;
-      
-      // edges of neighbour sides of B_FRONT corresponding to front's edges
-      EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
-      EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
-      for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i )
-        ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i]));
-      if ( ok )
-        ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP ));
 
+      // edges of adjacent sides of B_FRONT corresponding to front's edges
+      EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
+      EQuadEdge edgeOfAdj  [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
+      // first find all sides detectable w/o advanced analysis,
+      // then repeat the search, which then may pass without advanced analysis
+      for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
+      {
+        for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
+          if ( !block._side[i] ) // try to find 4 sides adjacent to front side
+            ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i],edgeOfAdj[i],advAnalys));
+        if ( ok || !advAnalys)
+          if ( !block._side[B_BACK] && block._side[B_TOP] ) // try to find back side by top one
+            ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP, advAnalys ));
+        if ( !advAnalys ) ok = true;
+      }
+      ok = block.isValid();
       if ( ok )
       {
         // check if just found block is same as one of previously found
         bool isSame = false;
         for ( int i = 1; i < _blocks.size() && !isSame; ++i )
-        {
-          _Block& prevBlock = _blocks[i-1];
-          isSame = true;
-          for ( int j = 0; j < 6 && isSame; ++j )
-            isSame = prevBlock.hasSide( block._side[ j ]);
-        }
+          isSame = ( block._corners == _blocks[i-1]._corners );
         ok = !isSame;
       }
 
       // count the found sides
-      _DUMP_(endl);
-      for (int i = 0; i < B_UNDEFINED; ++i )
+      _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
+      for (int i = 0; i < NB_BLOCK_SIDES; ++i )
       {
-        _DUMP_("** Block side "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
+        _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
         if ( block._side[ i ] )
         {
           if ( ok && i != B_FRONT)
@@ -497,22 +555,24 @@ namespace
             setSideBoundToBlock( *block._side[ i ]._side );
             nbBlockSides--;
           }
-          _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0));
-          _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0));    
-          _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1));
-          _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1));
+          _DUMP_("\t corners "<<
+                 block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
+                 block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
+                 block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
+                 block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
         }
         else
         {
-          _DUMP_("Not found"<<endl);
+          _DUMP_("\t not found"<<endl);
         }
       }
       if ( !ok )
-        block._side[0] = block._side[1] = block._side[2] =
-          block._side[3] = block._side[4] = block._side[5] = 0;
+        block.clear();
       else
         nbBlocks++;
     }
+    _DUMP_("Nb found blocks "<< nbBlocks <<endl);
+
     if ( nbBlocks == 0 && _error.empty() )
       return BAD_MESH_ERR;
 
@@ -525,7 +585,9 @@ namespace
    */
   //================================================================================
 
-  bool _Skin::fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad)
+  bool _Skin::fillSide( _BlockSide&             side,
+                        const SMDS_MeshElement* cornerQuad,
+                        const SMDS_MeshNode*    nCorner)
   {
     // Find out size of block side mesured in nodes and by the way find two rows
     // of nodes in two directions.
@@ -533,13 +595,6 @@ namespace
     int x, y, nbX, nbY;
     const SMDS_MeshElement* firstQuad = cornerQuad;
     {
-      // get corner node of cornerQuad
-      const SMDS_MeshNode* nCorner = 0;
-      for ( int i = 0; i < 4 && !nCorner; ++i )
-        if ( isCornerNode( firstQuad->GetNode(i)))
-          nCorner = firstQuad->GetNode(i);
-      if ( !nCorner ) return false;
-
       // get a node on block edge
       int iCorner = firstQuad->GetNodeIndex( nCorner );
       const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
@@ -611,7 +666,14 @@ namespace
       }
     }
 
-    return true;
+    // check side validity
+    bool ok =
+      side.getCornerFace( side.getCornerNode( 0, 0 )) &&
+      side.getCornerFace( side.getCornerNode( 1, 0 )) &&
+      side.getCornerFace( side.getCornerNode( 0, 1 )) &&
+      side.getCornerFace( side.getCornerNode( 1, 1 ));
+
+    return ok;
   }
 
   //================================================================================
@@ -622,7 +684,8 @@ namespace
 
   _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
                                            EQuadEdge sharedSideEdge1,
-                                           EQuadEdge sharedSideEdge2)
+                                           EQuadEdge sharedSideEdge2,
+                                           bool      withGeometricAnalysis)
   {
     _Block& block = _blocks.back();
     _OrientedBlockSide& side1 = block._side[ startBlockSide ];
@@ -637,67 +700,93 @@ namespace
     set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
 
     // exclude loaded sides of block from sidesOnEdge
-    int nbLoadedSides = 0;
-    for (int i = 0; i < B_UNDEFINED; ++i )
-    {
+    for (int i = 0; i < NB_BLOCK_SIDES; ++i )
       if ( block._side[ i ] )
-      {
-        nbLoadedSides++;
         sidesOnEdge.erase( block._side[ i ]._side );
-      }
-    }
+
     int nbSidesOnEdge = sidesOnEdge.size();
     _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
     if ( nbSidesOnEdge == 0 )
       return 0;
 
     _BlockSide* foundSide = 0;
-    if ( nbSidesOnEdge == 1 /*|| nbSidesOnEdge == 2 && nbLoadedSides == 1 */)
+    if ( nbSidesOnEdge == 1 )
     {
       foundSide = *sidesOnEdge.begin();
     }
     else
     {
-      // Select one of found sides most close to startBlockSide
-
-      // gravity center of already loaded block sides
-      gp_XYZ gc(0,0,0);
-      for (int i = 0; i < B_UNDEFINED; ++i )
-        if ( block._side[ i ] )
-          gc += block._side[ i ]._side->getGC();
-      gc /= nbLoadedSides;
-
-      gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
-      gp_Vec p1p2( p1, p2 );
-
-      const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
-      gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
-      gp_Vec side1Dir( p1, p1Op );
-      gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
-      _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
-             << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
-      
-      map < double , _BlockSide* > angleOfSide;
       set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
-      for (; sideIt != sidesOnEdge.end(); ++sideIt )
+      int nbLoadedSides = block.nbSides();
+      if ( nbLoadedSides > 1 )
       {
-        _BlockSide* sideI = *sideIt;
-        const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
-        gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
-        gp_Vec sideIDir( p1, p1Op );
-        // compute angle of (sideIDir projection to pln) and (X dir of pln)
-        gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
-        double angle = sideIDirProj.Angle( gp::DX2d() );
-        if ( angle < 0 ) angle += 2 * PI;
-        angleOfSide.insert( make_pair( angle, sideI ));
-        _DUMP_("  "<< sideI << " - side dir ("
-               << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
-               << " angle " << angle);
+        // Find the side having more than 2 corners common with already loaded sides
+        for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
+        {
+          _BlockSide* sideI = *sideIt;
+          int nbCommonCorners =
+            block._corners.count( sideI->getCornerNode(0,0)) +
+            block._corners.count( sideI->getCornerNode(1,0)) +
+            block._corners.count( sideI->getCornerNode(0,1)) +
+            block._corners.count( sideI->getCornerNode(1,1));
+          if ( nbCommonCorners > 2 )
+            foundSide = sideI;
+        }
+      }
+
+      if ( !foundSide )
+      {
+        if ( !withGeometricAnalysis ) return 0;
+
+        // Select one of found sides most close to startBlockSide
+
+        gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
+        gp_Vec p1p2( p1, p2 );
+
+        const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
+        gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
+        gp_Vec side1Dir( p1, p1Op );
+        gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
+        _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
+               << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
+
+        map < double , _BlockSide* > angleOfSide;
+        for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+        {
+          _BlockSide* sideI = *sideIt;
+          const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
+          gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
+          gp_Vec sideIDir( p1, p1Op );
+          // compute angle of (sideIDir projection to pln) and (X dir of pln)
+          gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
+          double angle = sideIDirProj.Angle( gp::DX2d() );
+          if ( angle < 0 ) angle += 2 * PI; // angle [0-2*PI]
+          angleOfSide.insert( make_pair( angle, sideI ));
+          _DUMP_("  "<< sideI << " - side dir ("
+                 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
+                 << " angle " << angle);
+        }
+        if ( nbLoadedSides == 1 )
+        {
+          double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first;
+          if ( angF > PI ) angF = 2*PI - angF;
+          if ( angL > PI ) angL = 2*PI - angL;
+          foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second;
+        }
+        else
+        {
+          gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
+          for (int i = 0; i < NB_BLOCK_SIDES; ++i )
+            if ( block._side[ i ] )
+              gc += block._side[ i ]._side->getGC();
+          gc /= nbLoadedSides;
+
+          gp_Vec gcDir( p1, gc );
+          gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
+          double gcAngle = gcDirProj.Angle( gp::DX2d() );
+          foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
+        }
       }
-      gp_Vec gcDir( p1, gc );
-      gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
-      double gcAngle = gcDirProj.Angle( gp::DX2d() );
-      foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
       _DUMP_("  selected "<< foundSide );
     }
 
@@ -759,9 +848,14 @@ namespace
       row2.push_back( n1 = oppositeNode( quad, i1 ));
     }
 
+    _NodeDescriptor nDesc( row1[1] );
+    if ( nDesc == _NodeDescriptor( row1[0] ))
+      return true; // row of 2 nodes
+
     // Find the rest nodes
     TIDSortedElemSet emptySet, avoidSet;
-    while ( !isCornerNode( n2 ))
+    //while ( !isCornerNode( n2 ))
+    while ( nDesc == _NodeDescriptor( n2 ))
     {
       avoidSet.clear(); avoidSet.insert( quad );
       quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
@@ -804,7 +898,31 @@ namespace
     return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
   }
 
-}
+  //================================================================================
+  /*!
+   * \brief Checks own validity
+   */
+  //================================================================================
+
+  bool _Block::isValid() const
+  {
+    bool ok = ( nbSides() == 6 );
+
+    // check only corners depending on side selection
+    EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
+    EQuadEdge edgeAdj [4] = { Q_TOP,    Q_RIGHT, Q_TOP, Q_RIGHT };
+    EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
+
+    for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
+    { 
+      SMESH_OrientedLink eBack = _side[ B_BACK      ].edge( edgeBack[i] );
+      SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
+      ok = ( eBack == eAdja );
+    }
+    return ok;
+  }
+
+} // namespace
 
 //=======================================================================
 //function : StdMeshers_HexaFromSkin_3D