Salome HOME
22360]: EDF SMESH: Body Fitting algorithm: incorporate edges
[modules/smesh.git] / src / StdMeshers / StdMeshers_Cartesian_3D.cxx
index d5435e7ee1b84197a567d3f98ff3902a8d37cf95..e880bb1a798f54fbb29fbe1c0aafe03729657b42 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -54,6 +54,8 @@
 #include <Geom2d_BSplineCurve.hxx>
 #include <Geom2d_BezierCurve.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+#include <GeomLib.hxx>
 #include <Geom_BSplineCurve.hxx>
 #include <Geom_BSplineSurface.hxx>
 #include <Geom_BezierCurve.hxx>
@@ -83,7 +85,7 @@
 #include <gp_Sphere.hxx>
 #include <gp_Torus.hxx>
 
-//#undef WITH_TBB
+#undef WITH_TBB
 #ifdef WITH_TBB
 #include <tbb/parallel_for.h>
 //#include <tbb/enumerable_thread_specific.h>
@@ -179,7 +181,7 @@ namespace
 
     B_IntersectPoint(): _node(NULL) {}
     void Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=0 ) const;
-    bool HasCommonFace( const B_IntersectPoint * other ) const;
+    int HasCommonFace( const B_IntersectPoint * other, int avoidFace=-1 ) const;
     bool IsOnFace( int faceID ) const;
     virtual ~B_IntersectPoint() {}
   };
@@ -277,7 +279,7 @@ namespace
   {
     vector< double >   _coords[3]; // coordinates of grid nodes
     gp_XYZ             _axes  [3]; // axis directions
-    vector< GridLine > _lines [3]; //  in 3 directions
+    vector< GridLine > _lines [3]; //    in 3 directions
     double             _tol, _minCellSize;
     gp_XYZ             _origin;
     gp_Mat             _invB; // inverted basis of _axes
@@ -289,6 +291,8 @@ namespace
     list< E_IntersectPoint >          _edgeIntP; // intersections with EDGEs
     TopTools_IndexedMapOfShape        _shapes;
 
+    SMESH_MesherHelper*               _helper;
+
     size_t CellIndex( size_t i, size_t j, size_t k ) const
     {
       return i + j*(_coords[0].size()-1) + k*(_coords[0].size()-1)*(_coords[1].size()-1);
@@ -442,8 +446,10 @@ namespace
     {
       const SMDS_MeshNode*       _node; // mesh node at hexahedron corner
       const B_IntersectPoint* _intPoint;
+      bool                _isUsedInFace;
 
-      _Node(const SMDS_MeshNode* n=0, const B_IntersectPoint* ip=0):_node(n), _intPoint(ip) {} 
+      _Node(const SMDS_MeshNode* n=0, const B_IntersectPoint* ip=0)
+        :_node(n), _intPoint(ip), _isUsedInFace(0) {} 
       const SMDS_MeshNode*    Node() const
       { return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
       const F_IntersectPoint* FaceIntPnt() const
@@ -463,9 +469,10 @@ namespace
           _intPoint->Add( ip->_faceIDs );
         }
       }
-      bool IsLinked( const B_IntersectPoint* other ) const
+      int IsLinked( const B_IntersectPoint* other,
+                    int                     avoidFace=-1 ) const // returns id of a common face
       {
-        return _intPoint && _intPoint->HasCommonFace( other );
+        return _intPoint ? _intPoint->HasCommonFace( other, avoidFace ) : 0;
       }
       bool IsOnFace( int faceID ) const // returns true if faceID is found
       {
@@ -502,9 +509,9 @@ namespace
         return _OrientedLink(&_link->_splits[_reverse ? NbResultLinks()-i-1 : i],_reverse);
       }
       _Node* FirstNode() const { return _link->_nodes[ _reverse ]; }
-      _Node* LastNode() const { return _link->_nodes[ !_reverse ]; }
+      _Node* LastNode()  const { return _link->_nodes[ !_reverse ]; }
       operator bool() const { return _link; }
-      vector< TGeomID > GetNotUsedFace(const set<TGeomID>& usedIDs ) const // returns supporting FACEs
+      vector< TGeomID > GetNotUsedFace(const set<TGeomID>& usedIDs ) const // returns supporting FACEs
       {
         vector< TGeomID > faces;
         const B_IntersectPoint *ip0, *ip1;
@@ -518,6 +525,11 @@ namespace
         }
         return faces;
       }
+      bool HasEdgeNodes() const
+      {
+        return ( dynamic_cast< const E_IntersectPoint* >( _link->_nodes[0]->_intPoint ) ||
+                 dynamic_cast< const E_IntersectPoint* >( _link->_nodes[1]->_intPoint ));
+      }
     };
     // --------------------------------------------------------------------------------
     struct _Face
@@ -593,13 +605,16 @@ namespace
                           vector< Hexahedron* >&  hexes,
                           int ijk[], int dIJK[] );
     bool findChain( _Node* n1, _Node* n2, _Face& quad, vector<_Node*>& chainNodes );
+    bool closePolygon( _Face* polygon, vector<_Node*>& chainNodes ) const;
     int  addElements(SMESH_MesherHelper& helper);
+    bool is1stNodeOut( _Link& link ) const;
     bool isInHole() const;
     bool checkPolyhedronSize() const;
     bool addHexa ();
     bool addTetra();
     bool addPenta();
     bool addPyra ();
+    bool debugDumpLink( _Link* link );
     _Node* FindEqualNode( vector< _Node >&        nodes,
                           const E_IntersectPoint* ip,
                           const double            tol2 )
@@ -725,7 +740,8 @@ namespace
         return ( ipBef->_transition == Trans_OUT );
       return ( ipBef->_transition != Trans_OUT );
     }
-    return prevIsOut; // _transition == Trans_TANGENT
+    // _transition == Trans_TANGENT
+    return !prevIsOut;
   }
   //================================================================================
   /*
@@ -749,15 +765,16 @@ namespace
   }
   //================================================================================
   /*
-   * Returns \c true if \a other B_IntersectPoint holds the same face ID
+   * Returns index of a common face if any, else zero
    */
-  bool B_IntersectPoint::HasCommonFace( const B_IntersectPoint * other ) const
+  int B_IntersectPoint::HasCommonFace( const B_IntersectPoint * other, int avoidFace ) const
   {
     if ( other )
       for ( size_t i = 0; i < other->_faceIDs.size(); ++i )
-        if ( IsOnFace( other->_faceIDs[i] ) )
-          return true;
-    return false;
+        if ( avoidFace != other->_faceIDs[i] &&
+             IsOnFace   ( other->_faceIDs[i] ))
+          return other->_faceIDs[i];
+    return 0;
   }
   //================================================================================
   /*
@@ -1582,11 +1599,11 @@ namespace
         _Link& link = _hexLinks[ iLink ];
         link._splits.clear();
         split._nodes[ 0 ] = link._nodes[0];
-        bool isOut = ( ! link._nodes[0]->Node() );
-        //int iEnd = link._intNodes.size() - bool( link._nodes[1]->_intPoint );
+        bool isOut = ( ! link._nodes[0]->Node() ); // is1stNodeOut( iLink );
+        bool checkTransition;
         for ( size_t i = 0; i < link._intNodes.size(); ++i )
         {
-          if ( link._intNodes[i].Node() )
+          if ( link._intNodes[i].Node() ) // intersection non-coinsident with a grid node
           {
             if ( split._nodes[ 0 ]->Node() && !isOut )
             {
@@ -1594,11 +1611,23 @@ namespace
               link._splits.push_back( split );
             }
             split._nodes[ 0 ] = &link._intNodes[i];
+            checkTransition = true;
+          }
+          else // FACE intersection coinsident with a grid node
+          {
+            checkTransition = ( link._nodes[0]->Node() );
           }
-          switch ( link._intNodes[i].FaceIntPnt()->_transition ) {
-          case Trans_OUT: isOut = true; break;
-          case Trans_IN : isOut = false; break;
-          default:; // isOut remains the same
+          if ( checkTransition )
+          {
+            switch ( link._intNodes[i].FaceIntPnt()->_transition ) {
+            case Trans_OUT: isOut = true; break;
+            case Trans_IN : isOut = false; break;
+            default:
+              if ( !link._intNodes[i].Node() && i == 0 )
+                isOut = is1stNodeOut( link );
+              else
+                  ; // isOut remains the same
+            }
           }
         }
         if ( link._nodes[ 1 ]->Node() && split._nodes[ 0 ]->Node() && !isOut )
@@ -1682,7 +1711,10 @@ namespace
           }
           break;
         }
-        default: // inside a hex
+        } // switch( nbFacets )
+
+        if ( nbFacets == 0 ||
+             _grid->_shapes( _edgeIntPnts[ iP ]->_shapeID ).ShapeType() == TopAbs_VERTEX )
         {
           equalNode = FindEqualNode( _vertexNodes, _edgeIntPnts[ iP ], tol2 );
           if ( equalNode ) {
@@ -1693,8 +1725,6 @@ namespace
             ++_nbIntNodes;
           }
         }
-        } // switch( nbFacets )
-
       } // loop on _edgeIntPnts
     }
     else if ( 3 < _nbCornerNodes && _nbCornerNodes < 8 ) // _nbIntNodes == 0
@@ -1744,9 +1774,10 @@ namespace
       return;
 
     _polygons.clear();
-    _polygons.reserve( 10 );
+    _polygons.reserve( 20 );
 
-    // create polygons from quadrangles and get their nodes
+    // Create polygons from quadrangles
+    // --------------------------------
 
     _Link polyLink;
     vector< _OrientedLink > splits;
@@ -1788,6 +1819,11 @@ namespace
             _vertexNodes.push_back( quad._edgeNodes[ iP ]);
         }
       }
+#ifdef _DEBUG_
+      for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
+        quad._edgeNodes[ iP ]._isUsedInFace = false;
+#endif
+      int nbUsedEdgeNodes = 0;
 
       while ( nbSplits > 0 )
       {
@@ -1815,8 +1851,8 @@ namespace
           n1 = split.FirstNode();
           if ( n1 != n2 )
           {
-            // try to connect to intersections with EDGES
-            if ( quad._edgeNodes.size() > 0  &&
+            // try to connect to intersections with EDGEs
+            if ( quad._edgeNodes.size() > nbUsedEdgeNodes  &&
                  findChain( n2, n1, quad, chainNodes ))
             {
               for ( size_t i = 1; i < chainNodes.size(); ++i )
@@ -1825,6 +1861,13 @@ namespace
                 polyLink._nodes[1] = chainNodes[i];
                 polygon->_polyLinks.push_back( polyLink );
                 polygon->_links.push_back( _OrientedLink( &polygon->_polyLinks.back() ));
+                nbUsedEdgeNodes += polyLink._nodes[1]->_isUsedInFace;
+              }
+              if ( chainNodes.back() != n1 )
+              {
+                n2 = chainNodes.back();
+                --iS;
+                continue;
               }
             }
             // try to connect to a split ending on the same FACE
@@ -1870,7 +1913,11 @@ namespace
 
         if ( nFirst != n2 ) // close a polygon
         {
-          findChain( n2, nFirst, quad, chainNodes );
+          if ( !findChain( n2, nFirst, quad, chainNodes ))
+          {
+            if ( !closePolygon( polygon, chainNodes ))
+              chainNodes.push_back( nFirst );
+          }
           for ( size_t i = 1; i < chainNodes.size(); ++i )
           {
             polyLink._nodes[0] = chainNodes[i-1];
@@ -1887,12 +1934,29 @@ namespace
         }
       } // while ( nbSplits > 0 )
 
+      if ( quad._edgeNodes.size() > nbUsedEdgeNodes )
+      {
+        // make _vertexNodes from not used _edgeNodes
+        const double tol = 0.05 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] );
+        for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
+        {
+          if ( quad._edgeNodes[ iP ]._isUsedInFace ) continue;
+          _Node* equalNode =
+            FindEqualNode( _vertexNodes, quad._edgeNodes[ iP ].EdgeIntPnt(), tol*tol );
+          if ( equalNode )
+            equalNode->Add( quad._edgeNodes[ iP ].EdgeIntPnt() );
+          else
+            _vertexNodes.push_back( quad._edgeNodes[ iP ]);
+        }
+      }
+
       if ( polygon->_links.size() < 3 )
         _polygons.pop_back();
 
     }  // loop on 6 sides of a hexahedron
 
-    // create polygons closing holes in a polyhedron
+    // Create polygons closing holes in a polyhedron
+    // ----------------------------------------------
 
     // add polygons to their links
     for ( size_t iP = 0; iP < _polygons.size(); ++iP )
@@ -1915,7 +1979,7 @@ namespace
           freeLinks.push_back( & polygon._links[ iL ]);
     }
     int nbFreeLinks = freeLinks.size();
-    if ( 0 < nbFreeLinks && nbFreeLinks < 3 ) return;
+    if ( nbFreeLinks > 0 && nbFreeLinks < 3 ) return;
 
     set<TGeomID> usedFaceIDs;
 
@@ -1929,7 +1993,8 @@ namespace
 
       _OrientedLink* curLink = 0;
       _Node*         curNode;
-      if ( !hasEdgeIntersections )
+      if (( !hasEdgeIntersections ) ||
+          ( nbFreeLinks < 4 && _vertexNodes.empty() ))
       {
         // get a remaining link to start from
         for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
@@ -1958,22 +2023,42 @@ namespace
         // get a remaining link to start from, one lying on minimal
         // nb of FACEs
         {
-          map< vector< TGeomID >, int > facesOfLink;
-          map< vector< TGeomID >, int >::iterator f2l;
+          vector< pair< TGeomID, int > > facesOfLink[3];
+          pair< TGeomID, int > faceOfLink( -1, -1 );
+          vector< TGeomID > faces;
           for ( size_t iL = 0; iL < freeLinks.size(); ++iL )
             if ( freeLinks[ iL ] )
             {
-              f2l = facesOfLink.insert
-                ( make_pair( freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs ), iL )).first;
-              if ( f2l->first.size() == 1 )
-                break;
+              faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs );
+              if ( faces.size() == 1 )
+              {
+                faceOfLink = make_pair( faces[0], iL );
+                if ( !freeLinks[ iL ]->HasEdgeNodes() )
+                  break;
+                facesOfLink[0].push_back( faceOfLink );
+              }
+              else if ( facesOfLink[0].empty() )
+              {
+                faceOfLink = make_pair(( faces.empty() ? -1 : faces[0]), iL );
+                facesOfLink[ 1 + faces.empty() ].push_back( faceOfLink );
+              }
             }
-          f2l = facesOfLink.begin();
-          if ( f2l->first.empty() )
-            return;
-          curFace = f2l->first[0];
-          curLink = freeLinks[ f2l->second ];
-          freeLinks[ f2l->second ] = 0;
+          for ( int i = 0; faceOfLink.second < 0 && i < 3; ++i )
+            if ( !facesOfLink[i].empty() )
+              faceOfLink = facesOfLink[i][0];
+
+          if ( faceOfLink.first < 0 ) // all faces used
+          {
+            for ( size_t i = 0; i < facesOfLink[2].size() && faceOfLink.first < 1; ++i )
+            {
+              curLink = freeLinks[ facesOfLink[2][i].second ];
+              faceOfLink.first = curLink->FirstNode()->IsLinked( curLink->FirstNode()->_intPoint );
+            }
+            usedFaceIDs.clear();
+          }
+          curFace = faceOfLink.first;
+          curLink = freeLinks[ faceOfLink.second ];
+          freeLinks[ faceOfLink.second ] = 0;
         }
         usedFaceIDs.insert( curFace );
         polygon._links.push_back( *curLink );
@@ -2043,12 +2128,15 @@ namespace
                 // TODO: to reorder _vertexNodes within polygon, if there are several ones
               }
           }
-          polyLink._nodes[0] = polygon._links[0].LastNode();
-          polyLink._nodes[1] = curNode;
-          polygon._polyLinks.push_back( polyLink );
-          polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() ));
-          freeLinks.push_back( &polygon._links.back() );
-          ++nbFreeLinks;
+          // if ( polygon._links.size() > 1 )
+          {
+            polyLink._nodes[0] = polygon._links[0].LastNode();
+            polyLink._nodes[1] = curNode;
+            polygon._polyLinks.push_back( polyLink );
+            polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() ));
+            freeLinks.push_back( &polygon._links.back() );
+            ++nbFreeLinks;
+          }
         }
 
       } // if there are intersections with EDGEs
@@ -2059,6 +2147,17 @@ namespace
 
       if ( polygon._links.size() == 2 )
       {
+        if ( freeLinks.back() == &polygon._links.back() )
+        {
+          freeLinks.back() = 0;
+          --nbFreeLinks;
+        }
+        vector< _Face*>& polygs1 = polygon._links.front()._link->_faces;
+        vector< _Face*>& polygs2 = polygon._links.back()._link->_faces;
+        _Face* polyg1 = ( polygs1.empty() ? 0 : polygs1[0] );
+        _Face* polyg2 = ( polygs2.empty() ? 0 : polygs2[0] );
+        if ( polyg1 ) polygs2.push_back( polyg1 );
+        if ( polyg2 ) polygs1.push_back( polyg2 );
         _polygons.pop_back();
       }
       else
@@ -2131,7 +2230,7 @@ namespace
         multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.begin();
         for ( ; ip != line._intPoints.end(); ++ip )
         {
-          //if ( !ip->_node ) continue;
+          // if ( !ip->_node ) continue; // intersection at a grid node
           lineInd.SetIndexOnLine( ip->_indexOnLine );
           for ( int iL = 0; iL < 4; ++iL ) // loop on 4 cells sharing a link
           {
@@ -2173,14 +2272,16 @@ namespace
       if ( hex )
       {
         intHexInd[ nbIntHex++ ] = i;
-        if ( hex->_nbIntNodes > 0 ) continue; // treat intersected hex later
+        if ( hex->_nbIntNodes > 0 || ! hex->_edgeIntPnts.empty())
+          continue; // treat intersected hex later
         this->init( hex->_i, hex->_j, hex->_k );
       }
       else
       {    
         this->init( i );
       }
-      if ( _nbCornerNodes == 8 && ( _nbBndNodes < _nbCornerNodes || !isInHole() ))
+      if (( _nbCornerNodes == 8 ) &&
+          ( _nbBndNodes < _nbCornerNodes || !isInHole() ))
       {
         // order of _hexNodes is defined by enum SMESH_Block::TShapeID
         SMDS_MeshElement* el =
@@ -2325,9 +2426,11 @@ namespace
         if ( iDirZ == 0 )
         {
           ip._point   = p1;
+          ip._shapeID = _grid->_shapes.Add( v1 );
           _grid->_edgeIntP.push_back( ip );
           if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 ))
             _grid->_edgeIntP.pop_back();
+          ip._shapeID = edgeID;
         }
         for ( int iP = 2; iP <= discret.NbPoints(); ++iP )
         {
@@ -2370,6 +2473,7 @@ namespace
         // add the 2nd vertex point to a hexahedron
         if ( iDirZ == 0 )
         {
+          ip._shapeID = _grid->_shapes.Add( v2 );
           ip._point = p1;
           _grid->ComputeUVW( p1, ip._uvw );
           locateValue( ijk[iDirX], ip._uvw[iDirX], _grid->_coords[iDirX], dIJK[iDirX], tol );
@@ -2378,6 +2482,7 @@ namespace
           _grid->_edgeIntP.push_back( ip );
           if ( !addIntersection( _grid->_edgeIntP.back(), hexes, ijk, d000 ))
             _grid->_edgeIntP.pop_back();
+          ip._shapeID = edgeID;
         }
       } // loop on 3 grid directions
     } // loop on EDGEs
@@ -2557,24 +2662,163 @@ namespace
   {
     chn.clear();
     chn.push_back( n1 );
-    bool found = false;
+    for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
+      if ( !quad._edgeNodes[ iP ]._isUsedInFace &&
+           n1->IsLinked( quad._edgeNodes[ iP ]._intPoint ) &&
+           n2->IsLinked( quad._edgeNodes[ iP ]._intPoint ))
+      {
+        chn.push_back( & quad._edgeNodes[ iP ]);
+        chn.push_back( n2 );
+        quad._edgeNodes[ iP ]._isUsedInFace = true;
+        return true;
+      }
+    bool found;
     do
     {
       found = false;
       for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP )
-        if (( std::find( ++chn.begin(), chn.end(), & quad._edgeNodes[iP]) == chn.end()) &&
-            chn.back()->IsLinked( quad._edgeNodes[ iP ]._intPoint ))
+        if ( !quad._edgeNodes[ iP ]._isUsedInFace &&
+             chn.back()->IsLinked( quad._edgeNodes[ iP ]._intPoint ))
         {
           chn.push_back( & quad._edgeNodes[ iP ]);
-          found = true;
+          found = quad._edgeNodes[ iP ]._isUsedInFace = true;
           break;
         }
-    } while ( found && chn.back() != n2 );
+    } while ( found && ! chn.back()->IsLinked( n2->_intPoint ) );
 
-    if ( chn.back() != n2 )
+    if ( chn.back() != n2 && chn.back()->IsLinked( n2->_intPoint ))
       chn.push_back( n2 );
 
-    return chn.size() > 2;
+    return chn.size() > 1;
+  }
+  //================================================================================
+  /*!
+   * \brief Try to heal a polygon whose ends are not connected
+   */
+  bool Hexahedron::closePolygon( _Face* polygon, vector<_Node*>& chainNodes ) const
+  {
+    int i = -1, nbLinks = polygon->_links.size();
+    if ( nbLinks < 3 )
+      return false;
+    vector< _OrientedLink > newLinks;
+    // find a node lying on the same FACE as the last one
+    _Node*   node = polygon->_links.back().LastNode();
+    int avoidFace = node->IsLinked( polygon->_links.back().FirstNode()->_intPoint );
+    for ( i = nbLinks - 2; i >= 0; --i )
+      if ( node->IsLinked( polygon->_links[i].FirstNode()->_intPoint, avoidFace ))
+        break;
+    if ( i >= 0 )
+    {
+      for ( ; i < nbLinks; ++i )
+        newLinks.push_back( polygon->_links[i] );
+    }
+    else
+    {
+      // find a node lying on the same FACE as the first one
+      node      = polygon->_links[0].FirstNode();
+      avoidFace = node->IsLinked( polygon->_links[0].LastNode()->_intPoint );
+      for ( i = 1; i < nbLinks; ++i )
+        if ( node->IsLinked( polygon->_links[i].LastNode()->_intPoint, avoidFace ))
+          break;
+      if ( i < nbLinks )
+        for ( nbLinks = i + 1, i = 0; i < nbLinks; ++i )
+          newLinks.push_back( polygon->_links[i] );
+    }
+    if ( newLinks.size() > 1 )
+    {
+      polygon->_links.swap( newLinks );
+      chainNodes.clear();
+      chainNodes.push_back( polygon->_links.back().LastNode() );
+      chainNodes.push_back( polygon->_links[0].FirstNode() );
+      return true;
+    }
+    return false;
+  }
+  //================================================================================
+  /*!
+   * \brief Checks transition at the 1st node of a link
+   */
+  bool Hexahedron::is1stNodeOut( _Link& link /*int iLink*/ ) const
+  {
+    // new version is for the case: tangent transition at the 1st node
+    bool isOut = false;
+    if ( link._intNodes.size() > 1 )
+    {
+      // check transition at the next intersection
+      switch ( link._intNodes[1].FaceIntPnt()->_transition ) {
+      case Trans_OUT: return false;
+      case Trans_IN : return true;
+      default: ; // tangent transition
+      }
+    }
+    gp_Pnt p1 = link._nodes[0]->Point();
+    gp_Pnt p2 = link._nodes[1]->Point();
+    gp_Pnt testPnt = 0.8 * p1.XYZ() + 0.2 * p2.XYZ();
+
+    TGeomID          faceID = link._intNodes[0]._intPoint->_faceIDs[0];
+    const TopoDS_Face& face = TopoDS::Face( _grid->_shapes( faceID ));
+    TopLoc_Location loc;
+    GeomAPI_ProjectPointOnSurf& proj =
+      _grid->_helper->GetProjector( face, loc, 0.1*_grid->_tol );
+    testPnt.Transform( loc );
+    proj.Perform( testPnt );
+    if ( proj.IsDone() &&
+         proj.NbPoints() > 0 &&
+         proj.LowerDistance() > _grid->_tol )
+    {
+      Quantity_Parameter u,v;
+      proj.LowerDistanceParameters( u,v );
+      gp_Dir normal;
+      if ( GeomLib::NormEstim( BRep_Tool::Surface( face, loc ),
+                               gp_Pnt2d( u,v ),
+                               0.1*_grid->_tol,
+                               normal ) < 3 )
+      {
+        if ( face.Orientation() == TopAbs_REVERSED )
+          normal.Reverse();
+        gp_Vec v( proj.NearestPoint(), testPnt );
+        return v * normal > 0;
+      }
+    }
+    //     if ( !_hexLinks[ iLink ]._nodes[0]->Node() ) // no node
+    //       return true;
+    //     if ( !_hexLinks[ iLink ]._nodes[0]->_intPoint ) // no intersection with geometry
+    //       return false;
+    //     switch ( _hexLinks[ iLink ]._nodes[0]->FaceIntPnt()->_transition ) {
+    //     case Trans_OUT: return true;
+    //     case Trans_IN : return false;
+    //     default: ; // tangent transition
+    //     }
+
+//     // ijk of a GridLine corresponding to the link
+//     int   iDir = iLink / 4;
+//     int indSub = iLink % 4;
+//     LineIndexer li = _grid->GetLineIndexer( iDir );
+//     li.SetIJK( _i,_j,_k );
+//     size_t lineIndex[4] = { li.LineIndex  (),
+//                             li.LineIndex10(),
+//                             li.LineIndex01(),
+//                             li.LineIndex11() };
+//     GridLine& line = _grid->_lines[ iDir ][ lineIndex[ indSub ]];
+
+//     // analyze transition of previous ip
+//     bool isOut = true;
+//     multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.begin();
+//     for ( ; ip != line._intPoints.end(); ++ip )
+//     {
+//       if ( &(*ip) == _hexLinks[ iLink ]._nodes[0]->_intPoint )
+//         break;
+//       switch ( ip->_transition ) {
+//       case Trans_OUT: isOut = true;
+//       case Trans_IN : isOut = false;
+//       default:;
+//       }
+//     }
+// #ifdef _DEBUG_
+//     if ( ip == line._intPoints.end() )
+//       cout << "BUG: Wrong GridLine. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl;
+// #endif
+    return isOut;
   }
   //================================================================================
   /*!
@@ -2728,7 +2972,9 @@ namespace
 
       // find a top node above the base node
       _Link* link = _polygons[0]._links[iL]._link;
-      ASSERT( link->_faces.size() > 1 );
+      //ASSERT( link->_faces.size() > 1 );
+      if ( link->_faces.size() < 2 )
+        return debugDumpLink( link );
       // a quadrangle sharing <link> with _polygons[0]
       _Face* quad = link->_faces[ bool( link->_faces[0] == & _polygons[0] )];
       for ( int i = 0; i < 4; ++i )
@@ -2757,7 +3003,9 @@ namespace
     nodes[2] = _polygons[0]._links[2].FirstNode();
 
     _Link* link = _polygons[0]._links[0]._link;
-    ASSERT( link->_faces.size() > 1 );
+    //ASSERT( link->_faces.size() > 1 );
+    if ( link->_faces.size() < 2 )
+      return debugDumpLink( link );
 
     // a triangle sharing <link> with _polygons[0]
     _Face* tria = link->_faces[ bool( link->_faces[0] == & _polygons[0] )];
@@ -2795,7 +3043,9 @@ namespace
 
       // find a top node above the base node
       _Link* link = _polygons[ iTri ]._links[iL]._link;
-      ASSERT( link->_faces.size() > 1 );
+      //ASSERT( link->_faces.size() > 1 );
+      if ( link->_faces.size() < 2 )
+        return debugDumpLink( link );
       // a quadrangle sharing <link> with a base triangle
       _Face* quad = link->_faces[ bool( link->_faces[0] == & _polygons[ iTri ] )];
       if ( quad->_links.size() != 4 ) return false;
@@ -2835,6 +3085,8 @@ namespace
 
     _Link* link = _polygons[iQuad]._links[0]._link;
     ASSERT( link->_faces.size() > 1 );
+    if ( link->_faces.size() < 2 )
+      return debugDumpLink( link );
 
     // a triangle sharing <link> with a base quadrangle
     _Face* tria = link->_faces[ bool( link->_faces[0] == & _polygons[ iQuad ] )];
@@ -2849,6 +3101,20 @@ namespace
 
     return false;
   }
+  //================================================================================
+  /*!
+   * \brief Dump a link and return \c false
+   */
+  bool Hexahedron::debugDumpLink( Hexahedron::_Link* link )
+  {
+#ifdef _DEBUG_
+    gp_Pnt p1 = link->_nodes[0]->Point(), p2 = link->_nodes[1]->Point();
+    cout << "BUG: not shared link. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl
+         << "n1 (" << p1.X() << ", "<< p1.Y() << ", "<< p1.Z() << " )" << endl
+         << "n2 (" << p2.X() << ", "<< p2.Y() << ", "<< p2.Z() << " )" << endl;
+#endif
+    return false;
+  }
 
   //================================================================================
   /*!
@@ -2964,9 +3230,12 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
 
   _computeCanceled = false;
 
+  SMESH_MesherHelper helper( theMesh );
+
   try
   {
     Grid grid;
+    grid._helper = &helper;
 
     vector< TopoDS_Shape > faceVec;
     {
@@ -2992,12 +3261,16 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
       shapeBox.Add( facesItersectors[i].GetFaceBndBox() );
 
       if ( _hyp->GetToAddEdges() )
+      {
+        helper.SetSubShape( faceVec[i] );
         for ( eExp.Init( faceVec[i], TopAbs_EDGE ); eExp.More(); eExp.Next() )
         {
           const TopoDS_Edge& edge = TopoDS::Edge( eExp.Current() );
-          if ( !SMESH_Algo::isDegenerated( edge ))
+          if ( !SMESH_Algo::isDegenerated( edge ) &&
+               !helper.IsRealSeam( edge ))
             edge2faceIDsMap[ grid._shapes.Add( edge )].push_back( facesItersectors[i]._faceID );
         }
+      }
     }
 
     getExactBndBox( faceVec, _hyp->GetAxisDirs(), shapeBox );
@@ -3037,7 +3310,6 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
     for ( size_t i = 0; i < facesItersectors.size(); ++i )
       facesItersectors[i].StoreIntersections();
 
-    SMESH_MesherHelper helper( theMesh );
     TopExp_Explorer solidExp (theShape, TopAbs_SOLID);
     helper.SetSubShape( solidExp.Current() );
     helper.SetElementsOnShape( true );