Fix SALOME_TESTS/Grids/smesh/imps_09/K0

Fix Triangulator for the case of self-intersecting but valid polygon

+ some fixes for #16469

1) ElementsOnShape: fix too high octree of classifiers in case of large tolerance
2) SMESH_MeshEditor::SewFreeBorder() SIGSEGV on over-constrained elements
3) Project(): adjust radius to avoid checking too many elements
4) Protect SMESH_Gen_i::Compute() from CORBA error in case of a removed mesh
5) smeshBuilder.Mesh.FindCoincidentNodesOnPart() - fix for a case of ID list

src/Controls/SMESH_Controls.cxx

@@ -4242,6 +4242,7 @@ struct ElementsOnShape::Classifier
   TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
   const TopoDS_Shape& Shape() const  { return myShape; }
   const Bnd_B3d* GetBndBox() const   { return & myBox; }
+  double Tolerance() const           { return myTol; }
   bool IsChecked()                   { return myFlags & theIsCheckedFlag; }
   bool IsSetFlag( int flag ) const   { return myFlags & flag; }
   void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
@@ -4458,10 +4459,9 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
     myOctree = new OctreeClassifier( myWorkClassifiers );
   }
 
-  SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
+  for ( int i = 0, nb = elem->NbNodes(); i < nb  && (isSatisfy == myAllNodesFlag); ++i )
-  while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
   {
-    SMESH_TNodeXYZ aPnt( aNodeItr->next() );
+    SMESH_TNodeXYZ aPnt( elem->GetNode( i ));
     centerXYZ += aPnt;
 
     isNodeOut = true;
@@ -4816,7 +4816,8 @@ void ElementsOnShape::OctreeClassifier::buildChildrenData()
   for ( int i = 0; i < nbChildren(); i++ )
   {
     OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
-    child->myIsLeaf = ( child->myClassifiers.size() <= 5  );
+    child->myIsLeaf = ( child->myClassifiers.size() <= 5  ||
+                        child->maxSize() < child->myClassifiers[0]->Tolerance() );
   }
 }
 

src/SMESH/SMESH_MeshEditor.cxx

@@ -7588,35 +7588,37 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
   theNodes.push_back( theSecondNode );
 
   const SMDS_MeshNode *nIgnore = theFirstNode, *nStart = theSecondNode;
-  TIDSortedElemSet foundElems;
+  //TIDSortedElemSet foundElems;
   bool needTheLast = ( theLastNode != 0 );
 
+  vector<const SMDS_MeshNode*> nodes;
+  
   while ( nStart != theLastNode ) {
     if ( nStart == theFirstNode )
       return !needTheLast;
 
-    // find all free border faces sharing form nStart
+    // find all free border faces sharing nStart
 
     list< const SMDS_MeshElement* > curElemList;
     list< const SMDS_MeshNode* >    nStartList;
     SMDS_ElemIteratorPtr invElemIt = nStart->GetInverseElementIterator(SMDSAbs_Face);
     while ( invElemIt->more() ) {
       const SMDS_MeshElement* e = invElemIt->next();
-      if ( e == curElem || foundElems.insert( e ).second ) {
+      //if ( e == curElem || foundElems.insert( e ).second ) // e can encounter twice in border
+      {
         // get nodes
-        int iNode = 0, nbNodes = e->NbNodes();
-        vector<const SMDS_MeshNode*> nodes( nbNodes+1 );
         nodes.assign( SMDS_MeshElement::iterator( e->interlacedNodesIterator() ),
                       SMDS_MeshElement::iterator() );
         nodes.push_back( nodes[ 0 ]);
 
         // check 2 links
+        int iNode = 0, nbNodes = nodes.size() - 1;
         for ( iNode = 0; iNode < nbNodes; iNode++ )
-          if (((nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
+          if ((( nodes[ iNode ] == nStart && nodes[ iNode + 1] != nIgnore ) ||
-               (nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
+               ( nodes[ iNode + 1] == nStart && nodes[ iNode ] != nIgnore )) &&
-              ControlFreeBorder( &nodes[ iNode ], e->GetID() ))
+              ( ControlFreeBorder( &nodes[ iNode ], e->GetID() )))
           {
-            nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart ? 1 : 0 )]);
+            nStartList.push_back( nodes[ iNode + ( nodes[ iNode ] == nStart )]);
             curElemList.push_back( e );
           }
       }
@@ -7668,7 +7670,7 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
         else if ( !contNodes[0].empty() && !contNodes[1].empty() ) {
           // choice: clear a worse one
           int iLongest = ( contNodes[0].size() < contNodes[1].size() ? 1 : 0 );
-          int iWorse = ( needTheLast ? 1 - iLongest : iLongest );
+          int   iWorse = ( needTheLast ? 1 - iLongest : iLongest );
           contNodes[ iWorse ].clear();
           contFaces[ iWorse ].clear();
         }
@@ -7682,10 +7684,8 @@ bool SMESH_MeshEditor::FindFreeBorder (const SMDS_MeshNode*             theFirst
       theNodes.pop_back(); // remove nIgnore
       theNodes.pop_back(); // remove nStart
       theFaces.pop_back(); // remove curElem
-      list< const SMDS_MeshNode* >::iterator nIt = cNL->begin();
+      theNodes.splice( theNodes.end(), *cNL );
-      list< const SMDS_MeshElement* >::iterator fIt = cFL->begin();
+      theFaces.splice( theFaces.end(), *cFL );
-      for ( ; nIt != cNL->end(); nIt++ ) theNodes.push_back( *nIt );
-      for ( ; fIt != cFL->end(); fIt++ ) theFaces.push_back( *fIt );
       return true;
 
     } // several continuations found
@@ -8026,6 +8026,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
     nIt[0] = nSide[0].begin(); eIt[0] = eSide[0].begin();
     nIt[1] = nSide[1].begin(); eIt[1] = eSide[1].begin();
 
+    // element can be split while iterating on border if it has two edges in the border
+    std::map< const SMDS_MeshElement* , const SMDS_MeshElement* > elemReplaceMap;
+    std::map< const SMDS_MeshElement* , const SMDS_MeshElement* >::iterator elemReplaceMapIt;
+
     TElemOfNodeListMap insertMap;
     TElemOfNodeListMap::iterator insertMapIt;
     // insertMap is
@@ -8073,12 +8077,15 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
         const SMDS_MeshNode*    nIns = *nIt [ 1 - intoBord ];
         if ( intoBord == 1 ) {
           // move node of the border to be on a link of elem of the side
-          gp_XYZ p1 (n1->X(), n1->Y(), n1->Z());
+          SMESH_NodeXYZ p1( n1 ), p2( n2 );
-          gp_XYZ p2 (n2->X(), n2->Y(), n2->Z());
           double ratio = du / ( param[ 1 ][ i[1] ] - param[ 1 ][ i[1]-1 ]);
           gp_XYZ p = p2 * ( 1 - ratio ) + p1 * ratio;
           GetMeshDS()->MoveNode( nIns, p.X(), p.Y(), p.Z() );
         }
+        elemReplaceMapIt = elemReplaceMap.find( elem );
+        if ( elemReplaceMapIt != elemReplaceMap.end() )
+          elem = elemReplaceMapIt->second;
+
         insertMapIt = insertMap.find( elem );
         bool  notFound = ( insertMapIt == insertMap.end() );
         bool otherLink = ( !notFound && (*insertMapIt).second.front() != n1 );
@@ -8101,8 +8108,10 @@ SMESH_MeshEditor::SewFreeBorder (const SMDS_MeshNode* theBordFirstNode,
             UpdateVolumes(n12, n22, nodeList);
           }
           // 3. find an element appeared on n1 and n2 after the insertion
-          insertMap.erase( elem );
+          insertMap.erase( insertMapIt );
-          elem = findAdjacentFace( n1, n2, 0 );
+          const SMDS_MeshElement* elem2 = findAdjacentFace( n1, n2, 0 );
+          elemReplaceMap.insert( std::make_pair( elem, elem2 ));
+          elem = elem2;
         }
         if ( notFound || otherLink ) {
           // add element and nodes of the side into the insertMap

src/SMESHUtils/SMESH_MeshAlgos.cxx

@@ -1254,18 +1254,43 @@ gp_XYZ SMESH_ElementSearcherImpl::Project(const gp_Pnt&            point,
   gp_XYZ p = point.XYZ();
   ElementBndBoxTree* ebbLeaf = ebbTree->getLeafAtPoint( p );
   const Bnd_B3d* box = ebbLeaf ? ebbLeaf->getBox() : ebbTree->getBox();
-  double radius = ( box->CornerMax() - box->CornerMin() ).Modulus();
+  gp_XYZ pMin = box->CornerMin(), pMax = box->CornerMax();
+  double radius = Precision::Infinite();
+  if ( ebbLeaf || !box->IsOut( p ))
+  {
+    for ( int i = 1; i <= 3; ++i )
+    {
+      double d = 0.5 * ( pMax.Coord(i) - pMin.Coord(i) );
+      if ( d > Precision::Confusion() )
+        radius = Min( d, radius );
+    }
+    if ( !ebbLeaf )
+      radius /= ebbTree->getHeight( /*full=*/true );
+  }
+  else // p outside of box
+  {
+    for ( int i = 1; i <= 3; ++i )
+    {
+      double d = 0;
+      if ( point.Coord(i) < pMin.Coord(i) )
+        d = pMin.Coord(i) - point.Coord(i);
+      else if ( point.Coord(i) > pMax.Coord(i) )
+        d = point.Coord(i) - pMax.Coord(i);
+      if ( d > Precision::Confusion() )
+        radius = Min( d, radius );
+    }
+  }
 
   ElementBndBoxTree::TElemSeq elems;
   ebbTree->getElementsInSphere( p, radius, elems );
   while ( elems.empty() && radius < 1e100 )
   {
-    radius *= 1.5;
+    radius *= 1.1;
     ebbTree->getElementsInSphere( p, radius, elems );
   }
   gp_XYZ proj, bestProj;
   const SMDS_MeshElement* elem = 0;
-  double minDist = 2 * radius;
+  double minDist = Precision::Infinite();
   ElementBndBoxTree::TElemSeq::iterator e = elems.begin();
   for ( ; e != elems.end(); ++e )
   {
@@ -2347,28 +2372,16 @@ void SMESH_MeshAlgos::Get1DBranches( SMDS_ElemIteratorPtr theEdgeIt,
 
     if ( nbBranches == 2 && !startIsBranchEnd ) // join two branches starting at the same node
     {
-      if ( nodeBranches[0].back() == nodeBranches[1].back() )
+      std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() );
-      {
+      nodeBranches[0].pop_back();
-        // it is a closed branch, keep theStartNode first
+      nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
-        nodeBranches[0].pop_back();
+      nodeBranches[0].insert( nodeBranches[0].end(),
-        nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
+                              nodeBranches[1].begin(), nodeBranches[1].end() );
-        nodeBranches[0].insert( nodeBranches[0].end(),
+
-                                nodeBranches[1].rbegin(), nodeBranches[1].rend() );
+      std::reverse( branches[0].begin(), branches[0].end() );
-        branches[0].reserve( branches[0].size() + branches[1].size() );
+      branches[0].reserve( branches[0].size() + branches[1].size() );
-        branches[0].insert( branches[0].end(), branches[1].rbegin(), branches[1].rend() );
+      branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() );
-      }
-      else
-      {
-        std::reverse( nodeBranches[0].begin(), nodeBranches[0].end() );
-        nodeBranches[0].pop_back();
-        nodeBranches[0].reserve( nodeBranches[0].size() + nodeBranches[1].size() );
-        nodeBranches[0].insert( nodeBranches[0].end(),
-                                nodeBranches[1].begin(), nodeBranches[1].end() );
 
-        std::reverse( branches[0].begin(), branches[0].end() );
-        branches[0].reserve( branches[0].size() + branches[1].size() );
-        branches[0].insert( branches[0].end(), branches[1].begin(), branches[1].end() );
-      }
       nodeBranches[1].clear();
       branches[1].clear();
     }

src/SMESHUtils/SMESH_MeshAlgos.hxx

@@ -443,28 +443,12 @@ namespace SMESH_MeshAlgos
 
     bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes );
 
-    /*!
+    struct PolyVertex;
-     * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
-     */
-    struct PolyVertex
-    {
-      SMESH_NodeXYZ _nxyz;
-      size_t        _index;
-      gp_XY         _xy;
-      PolyVertex*   _prev;
-      PolyVertex*   _next;
-
-      void   SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index );
-      void   GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const;
-      double TriaArea() const;
-      bool   IsInsideTria( const PolyVertex* v );
-      PolyVertex* Delete();
-    };
     struct Optimizer;
+    struct Data;
 
-    std::vector< PolyVertex > _pv;
+    Data*      _data;
-    std::vector< size_t >     _nodeIndex;
+    Optimizer* _optimizer;
-    Optimizer*                _optimizer;
   };
 
   // structure used in MakePolyLine() to define a cutting plane

src/SMESHUtils/SMESH_Offset.cxx

@@ -141,11 +141,11 @@ namespace
     EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector<const SMDS_MeshNode *>() ) {}
     void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; }
     bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; }
-    bool Contains( const SMDS_MeshNode* n ) const
+    size_t Contains( const SMDS_MeshNode* n ) const
     {
       for ( size_t i = 0; i < myLinks.size(); ++i )
-        if ( myLinks[i]->myNode1 == n ) return true;
+        if ( myLinks[i]->myNode1 == n ) return i + 1;
-      return false;
+      return 0;
     }
     virtual int NbNodes() const { return myLinks.size(); }
     virtual SMDS_ElemIteratorPtr nodesIterator() const
@@ -224,6 +224,24 @@ namespace
         myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0;
       loop->Clear();
     }
+    void Join( EdgeLoop& loop1, size_t iAfterConcact,
+               EdgeLoop& loop2, size_t iFromEdge2 )
+    {
+      std::vector< const EdgePart* > linksAfterContact( loop1.myLinks.begin() + iAfterConcact,
+                                                        loop1.myLinks.end() );
+      loop1.myLinks.reserve( loop2.myLinks.size() + loop1.myLinks.size() );
+      loop1.myLinks.resize( iAfterConcact );
+      loop1.myLinks.insert( loop1.myLinks.end(),
+                            loop2.myLinks.begin() + iFromEdge2, loop2.myLinks.end() );
+      loop1.myLinks.insert( loop1.myLinks.end(),
+                            loop2.myLinks.begin(), loop2.myLinks.begin() + iFromEdge2 );
+      loop1.myLinks.insert( loop1.myLinks.end(),
+                            linksAfterContact.begin(), linksAfterContact.end() );
+      loop1.myIsBndConnected = loop2.myIsBndConnected;
+      loop2.Clear();
+      for ( size_t iE = 0; iE < loop1.myLinks.size(); ++iE )
+        myLoopOfEdge[ Index( *loop1.myLinks[ iE ] )] = & loop1;
+    }
     size_t    Index( const EdgePart& edge ) const { return &edge - myEdge0; }
     EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; }
   };
@@ -2558,7 +2576,8 @@ namespace
   //================================================================================
   /*!
    * \brief Remove loops that are not connected to boundary edges of myFace by
-   *        adding edges connecting these loops to the boundary
+   *        adding edges connecting these loops to the boundary.
+   *        Such loops must be removed as they form polygons with ring topology.
    */
   //================================================================================
 
@@ -2607,14 +2626,49 @@ namespace
     while ( prevNbReached < nbReachedLoops );
 
 
-    // add links connecting internal loops with the boundary ones
 
     for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL )
     {
       EdgeLoop& loop = theLoops.myLoops[ iL ];
-      if ( loop.myIsBndConnected )
+      if ( loop.myIsBndConnected || loop.myLinks.size() == 0 )
         continue;
 
+      if ( loop.myHasPending )
+      {
+        // try to join the loop to another one, with which it contacts at a node
+
+        // look for a node where the loop reverses
+        const EdgePart* edgePrev = loop.myLinks.back();
+        for ( size_t iE = 0; iE < loop.myLinks.size(); edgePrev = loop.myLinks[ iE++ ] )
+        {
+          if ( !edgePrev->IsTwin( *loop.myLinks[ iE ]))
+            continue;
+          const SMDS_MeshNode* reverseNode = edgePrev->myNode2;
+
+          // look for a loop including reverseNode
+          size_t iContactEdge2; // index(+1) of edge starting at reverseNode
+          for ( size_t iL2 = 0; iL2 < theLoops.myNbLoops; ++iL2 )
+          {
+            if ( iL == iL2 )
+              continue;
+            EdgeLoop& loop2 = theLoops.myLoops[ iL2 ];
+            if ( ! ( iContactEdge2 = loop2.Contains( reverseNode )))
+              continue;
+
+            // insert loop2 into the loop
+            theLoops.Join( loop, iE, loop2, iContactEdge2 - 1 );
+            break;
+          }
+          if ( loop.myIsBndConnected )
+            break;
+        }
+
+        if ( loop.myIsBndConnected )
+          continue;
+      }
+
+      // add links connecting internal loops with the boundary ones
+
       // find a pair of closest nodes
       const SMDS_MeshNode *closestNode1, *closestNode2;
       double minDist = 1e100;
@@ -2689,7 +2743,7 @@ namespace
 
     while ( !theLoops.AllEdgesUsed() )
     {
-      theLoops.AddNewLoop();
+      EdgeLoop& loop = theLoops.AddNewLoop();
 
       // add 1st edge to a new loop
       size_t i1;
@@ -2720,7 +2774,7 @@ namespace
         // choose among candidates
         if ( theLoops.myCandidates.size() == 0 )
         {
-          theLoops.GetLoopOf( lastEdge )->myHasPending = true;
+          loop.myHasPending = bool( twinEdge );
           lastEdge = twinEdge;
         }
         else if ( theLoops.myCandidates.size() == 1 )
@@ -2751,6 +2805,10 @@ namespace
       }
       while ( lastNode != firstNode );
 
+
+      if ( twinEdge == & myLinks[ i1 ])
+        loop.myHasPending = true;
+
     } // while ( !theLoops.AllEdgesUsed() )
 
     return;

src/SMESHUtils/SMESH_PolyLine.cxx

@@ -52,7 +52,17 @@ namespace
     int                     mySrcPntInd; //!< start point index
     TIDSortedElemSet        myElemSet, myAvoidSet;
 
-    Path(): myLength(0.0), myFace(0) {}
+    Path(const SMDS_MeshElement* face=0, int srcInd=-1):
+      myLength(0.0), myFace(face), mySrcPntInd( srcInd ) {}
+
+    void CopyNoPoints( const Path& other );
+
+    bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths = 0 );
+
+    bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
+                         const gp_XYZ&           plnNorm,
+                         const gp_XYZ&           plnOrig,
+                         std::vector< Path >*    paths);
 
     bool SetCutAtCorner( const SMESH_NodeXYZ&    cornerNode,
                          const SMDS_MeshElement* face,
@@ -61,8 +71,6 @@ namespace
 
     void AddPoint( const gp_XYZ& p );
 
-    bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig );
-
     bool ReachSamePoint( const Path& other );
 
     static void Remove( std::vector< Path > & paths, size_t& i );
@@ -80,6 +88,25 @@ namespace
              myFace == other.myFace );
   }
 
+  //================================================================================
+  /*!
+   * \brief Copy data except points
+   */
+  //================================================================================
+
+  void Path::CopyNoPoints( const Path& other )
+  {
+    myLength    = other.myLength;
+    mySrcPntInd = other.mySrcPntInd;
+    myFace      = other.myFace;
+    myNode1     = other.myNode1;
+    myNode2     = other.myNode2;
+    myNodeInd1  = other.myNodeInd1;
+    myNodeInd2  = other.myNodeInd2;
+    myDot1      = other.myDot1;
+    myDot2      = other.myDot2;
+  }
+
   //================================================================================
   /*!
    * \brief Remove a path from a vector
@@ -93,16 +120,8 @@ namespace
       size_t j = paths.size() - 1; // last item to be removed
       if ( i < j )
       {
+        paths[ i ].CopyNoPoints ( paths[ j ]);
         paths[ i ].myPoints.swap( paths[ j ].myPoints );
-        paths[ i ].myLength    = paths[ j ].myLength;
-        paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd;
-        paths[ i ].myFace      = paths[ j ].myFace;
-        paths[ i ].myNode1     = paths[ j ].myNode1;
-        paths[ i ].myNode2     = paths[ j ].myNode2;
-        paths[ i ].myNodeInd1  = paths[ j ].myNodeInd1;
-        paths[ i ].myNodeInd2  = paths[ j ].myNodeInd2;
-        paths[ i ].myDot1      = paths[ j ].myDot1;
-        paths[ i ].myDot2      = paths[ j ].myDot2;
       }
     }
     paths.pop_back();
@@ -110,6 +129,62 @@ namespace
       --i;
   }
 
+  //================================================================================
+  /*!
+   * \brief Try to extend self by a point located at a node.
+   *        Return a success flag.
+   */
+  //================================================================================
+
+  bool Path::SetCutAtCorner( const SMESH_NodeXYZ&  cornerNode,
+                             const gp_XYZ&         plnNorm,
+                             const gp_XYZ&         plnOrig,
+                             std::vector< Path > * paths )
+  {
+    bool ok = false;
+    const bool isContinuation = myFace; // extend this path or find all possible paths?
+    const SMDS_MeshElement* lastFace = myFace;
+
+    myAvoidSet.clear();
+
+    SMDS_ElemIteratorPtr fIt = cornerNode->GetInverseElementIterator(SMDSAbs_Face);
+    while ( fIt->more() )
+    {
+      Path path( lastFace, mySrcPntInd );
+      if ( !path.SetCutAtCorner( cornerNode, fIt->next(), plnNorm, plnOrig ))
+        continue;
+
+      if ( !myAvoidSet.insert( path.myNode1.Node() ).second ||
+           !myAvoidSet.insert( path.myNode2.Node() ).second )
+        continue;
+
+      if ( isContinuation )
+      {
+        if ( ok ) // non-manifold continuation
+        {
+          path.myPoints = myPoints;
+          path.myLength = myLength;
+          path.AddPoint( cornerNode );
+          paths->push_back( path );
+        }
+        else
+        {
+          double len = myLength;
+          this->CopyNoPoints( path );
+          this->myLength = len;
+          this->AddPoint( path.myPoints.back() );
+        }
+      }
+      else
+      {
+        paths->push_back( path );
+      }
+      ok = true;
+    }
+
+    return ok;
+  }
+
   //================================================================================
   /*!
    * \brief Store a point that is at a node of a face if the face is intersected by plane.
@@ -169,11 +244,14 @@ namespace
    * \brief Try to find the next point
    *  \param [in] plnNorm - cutting plane normal
    *  \param [in] plnOrig - cutting plane origin
+   *  \param [in] paths - all paths
    */
   //================================================================================
 
-  bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig )
+  bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig, std::vector< Path > * paths )
   {
+    bool ok = false;
+
     int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
     if ( myNodeInd2 == nodeInd3 )
       nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
@@ -195,20 +273,13 @@ namespace
     }
     else if ( dot3 == 0. )
     {
-      SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face);
+      ok = SetCutAtCorner( node3, plnNorm, plnOrig, paths );
-      while ( fIt->more() )
+      return ok;
-        if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig ))
-          return true;
-      return false;
     }
     else if ( myDot2 == 0. )
     {
-      SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner()
+      ok = SetCutAtCorner( myNode2, plnNorm, plnOrig, paths );
-      SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face);
+      return ok;
-      while ( fIt->more() )
-        if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig ))
-          return true;
-      return false;
     }
 
     double r = Abs( myDot1 / ( myDot2 - myDot1 ));
@@ -216,10 +287,32 @@ namespace
 
     myAvoidSet.clear();
     myAvoidSet.insert( myFace );
-    myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
+    const SMDS_MeshElement* nextFace;
-                                             myElemSet,   myAvoidSet,
+    int ind1, ind2;
-                                             &myNodeInd1, &myNodeInd2 );
+    while (( nextFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
-    return myFace;
+                                                        myElemSet, myAvoidSet,
+                                                        &ind1, &ind2 )))
+    {
+      if ( ok ) // non-manifold continuation
+      {
+        paths->push_back( *this );
+        paths->back().myFace = nextFace;
+        paths->back().myNodeInd1 = ind1;
+        paths->back().myNodeInd2 = ind2;
+      }
+      else
+      {
+        myFace = nextFace;
+        myNodeInd1 = ind1;
+        myNodeInd2 = ind2;
+      }
+      ok = true;
+      if ( !paths )
+        break;
+      myAvoidSet.insert( nextFace );
+    }
+
+    return ok;
   }
 
   //================================================================================
@@ -263,6 +356,13 @@ namespace
     {
       SMESH_MeshAlgos::PolySegment& polySeg = mySegments[ iSeg ];
 
+      if ( ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ).SquareModulus() == 0 )
+      {
+        myPaths[ iSeg ].AddPoint( polySeg.myXYZ[0] );
+        myPaths[ iSeg ].AddPoint( polySeg.myXYZ[1] );
+        return;
+      }
+
       // the cutting plane
       gp_XYZ plnNorm = ( polySeg.myXYZ[0] - polySeg.myXYZ[1] ) ^ polySeg.myVector.XYZ();
       gp_XYZ plnOrig = polySeg.myXYZ[1];
@@ -275,14 +375,13 @@ namespace
 
       for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
       {
-        Path path;
+        Path path( 0, iP );
-        path.mySrcPntInd = iP;
         size_t nbPaths = paths.size();
 
         if ( polySeg.myFace[ iP ]) // the end point lies on polySeg.myFace[ iP ]
         {
           // check coincidence of polySeg.myXYZ[ iP ] with nodes
-          const double tol = 1e-20;
+          const double tol = 1e-17;
           SMESH_NodeXYZ nodes[4];
           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
           {
@@ -292,6 +391,11 @@ namespace
           }
           nodes[ 3 ] = nodes[ 0 ];
 
+          double dot[ 4 ];
+          for ( int i = 0; i < 3; ++i )
+            dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
+          dot[ 3 ] = dot[ 0 ];
+
           // check coincidence of polySeg.myXYZ[ iP ] with edges
           for ( int i = 0; i < 3 &&  !polySeg.myNode1[ iP ]; ++i )
           {
@@ -300,16 +404,35 @@ namespace
             {
               polySeg.myNode1[ iP ] = nodes[ i     ].Node();
               polySeg.myNode2[ iP ] = nodes[ i + 1 ].Node();
+
+              int i3 = ( i + 2 ) % 3;
+              if ( dot[ i   ] * dot [ i3 ] > 0 &&
+                   dot[ i+1 ] * dot [ i3 ] > 0 ) // point iP is inside a neighbor triangle
+              {
+                path.myAvoidSet.insert( polySeg.myFace[ iP ]);
+                const SMDS_MeshElement* face2 = 
+                  SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
+                                                  polySeg.myNode2[ iP ],
+                                                  path.myElemSet,
+                                                  path.myAvoidSet );
+                if ( face2 )
+                  polySeg.myFace[ iP ] = face2;
+                else
+                  ;// ??
+                for ( int i = 0; i < 3; ++i )
+                {
+                  nodes[ i ] = polySeg.myFace[ iP ]->GetNode( i );
+                  dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
+                }
+                dot[ 3 ] = dot[ 0 ];
+                polySeg.myNode1[ iP ] = polySeg.myNode2[ iP ] = 0;
+                break;
+              }
             }
           }
 
           if ( !polySeg.myNode1[ iP ] ) // polySeg.myXYZ[ iP ] is within polySeg.myFace[ iP ]
           {
-            double dot[ 4 ];
-            for ( int i = 0; i < 3; ++i )
-              dot[ i ] = plnNorm * ( nodes[ i ] - plnOrig );
-            dot[ 3 ] = dot[ 0 ];
-
             int iCut = 0; // index of a cut edge
             if      ( dot[ 1 ] * dot[ 2 ] < 0. ) iCut = 1;
             else if ( dot[ 2 ] * dot[ 3 ] < 0. ) iCut = 2;
@@ -364,27 +487,45 @@ namespace
             path.AddPoint( polySeg.myXYZ[ iP ]);
             path.myAvoidSet.insert( path.myFace );
             paths.push_back( path );
+            std::swap( polySeg.myNode1[ iP ], polySeg.myNode2[ iP ]);
           }
           if ( nbPaths == paths.size() )
             throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
                                      << " in a PolySegment " << iSeg );
-        }
+
-        else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
+          if ( path.myDot1 == 0. &&
-        {
+               path.myDot2 == 0. &&
-          std::set<const SMDS_MeshNode* > nodes;
+               paths.size() - nbPaths >= 2 ) // use a face non-parallel to the plane
-          SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
-          while ( fIt->more() )
           {
-            path.myPoints.clear();
+            const SMDS_MeshElement* goodFace = 0;
-            if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig ))
+            for ( size_t j = nbPaths; j < paths.size(); ++j )
             {
-              if (( path.myDot1 * path.myDot2 != 0 ) ||
+              path = paths[j];
-                  ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second ))
+              if ( path.Extend( plnNorm, plnOrig ))
-                paths.push_back( path );
+                goodFace = paths[j].myFace;
+              else
+                paths[j].myFace = 0;
             }
+            if ( !goodFace )
+              throw SALOME_Exception ( SMESH_Comment("Cant move from point ") << iP+1
+                                       << " of a PolySegment " << iSeg );
+            for ( size_t j = nbPaths; j < paths.size(); ++j )
+              if ( !paths[j].myFace )
+              {
+                paths[j].myFace = goodFace;
+                paths[j].myNodeInd1 = goodFace->GetNodeIndex( paths[j].myNode1.Node() );
+                paths[j].myNodeInd2 = goodFace->GetNodeIndex( paths[j].myNode2.Node() );
+              }
           }
         }
 
+        else if ( polySeg.myNode1[ iP ] ) // the end point is at a node
+        {
+          path.myFace = 0;
+          path.SetCutAtCorner( polySeg.myNode1[ iP ], plnNorm, plnOrig, &paths );
+        }
+
+
         // look for a one-segment path
         for ( size_t i = 0; i < nbPaths; ++i )
           for ( size_t j = nbPaths; j < paths.size(); ++j )
@@ -394,7 +535,9 @@ namespace
               myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
               paths.clear();
             }
-      }
+
+      }  // loop on the polySeg end points to initialize all possible paths
+
 
       // 2) extend paths and compose the shortest one connecting the two points
 
@@ -405,8 +548,8 @@ namespace
         for ( size_t i = 0; i < paths.size(); ++i )
         {
           Path& path = paths[ i ];
-          if ( !path.Extend( plnNorm, plnOrig ) ||        // path reached a mesh boundary
+          if ( !path.Extend( plnNorm, plnOrig, &paths ) || // path reached a mesh boundary
-               path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
+               path.myLength > myPaths[ iSeg ].myLength )  // path is longer than others
           {
             Path::Remove( paths, i );
             continue;
@@ -428,8 +571,10 @@ namespace
                                   paths[j].myPoints.rbegin(),
                                   paths[j].myPoints.rend() );
               }
+              if ( i < j ) std::swap( i, j );
               Path::Remove( paths, i );
               Path::Remove( paths, j );
+              break;
             }
           }
         }
@@ -504,7 +649,7 @@ void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh*                            theMes
     gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
 
     isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
-    if ( !isVectorOK[ iSeg ])
+    if ( !isVectorOK[ iSeg ] && ( p1 - p2 ).SquareModulus() > 0. )
     {
       gp_XYZ pMid = 0.5 * ( p1 + p2 );
       const SMDS_MeshElement* face;
@@ -512,14 +657,35 @@ void SMESH_MeshAlgos::MakePolyLine( SMDS_Mesh*                            theMes
       polySeg.myVector       = polySeg.myMidProjPoint.XYZ() - pMid;
 
       gp_XYZ faceNorm;
-      SMESH_MeshAlgos::FaceNormal( face, faceNorm );
+      SMESH_MeshAlgos::FaceNormal( face, faceNorm, /*normalized=*/false );
 
-      if ( polySeg.myVector.Magnitude() < Precision::Confusion() ||
+      const double tol = Precision::Confusion();
-           polySeg.myVector * faceNorm  < Precision::Confusion() )
+      if ( polySeg.myVector.Magnitude() < tol || polySeg.myVector * faceNorm  < tol )
       {
         polySeg.myVector = faceNorm;
         polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
       }
+      plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+      if ( plnNorm.SquareModulus() == 0 ) // p1-p2 perpendicular to mesh
+      {
+        double radius = faceNorm.Modulus();
+        std::vector< const SMDS_MeshElement* > foundElems;
+        while ( plnNorm.SquareModulus() == 0  &&  radius < 1e200 )
+        {
+          foundElems.clear();
+          searcher->GetElementsInSphere( p1, radius, SMDSAbs_Face, foundElems );
+          searcher->GetElementsInSphere( p2, radius, SMDSAbs_Face, foundElems );
+          radius *= 2;
+          polySeg.myVector.SetCoord( 0,0,0 );
+          for ( size_t i = 0; i < foundElems.size(); ++i )
+          {
+            SMESH_MeshAlgos::FaceNormal( foundElems[i], faceNorm );
+            polySeg.myVector += faceNorm / foundElems.size();
+          }
+          plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+        }
+      }
+      
     }
     else
     {

src/SMESHUtils/SMESH_Triangulate.cxx

@@ -60,13 +60,47 @@ namespace
     }
     bool operator<(const Node& other) const { return _triaIndex < other._triaIndex; }
   };
-  typedef boost::container::flat_set< Node > TNodeSet;
+  typedef boost::container::flat_set< Node >                 TriaNodeSet;
 
 }
+/*!
+ * \brief Vertex of a polygon. Together with 2 neighbor Vertices represents a triangle
+ */
+struct Triangulate::PolyVertex
+{
+  SMESH_NodeXYZ _nxyz;
+  size_t        _index;
+  gp_XY         _xy;
+  PolyVertex*   _prev;
+  PolyVertex*   _next;
+
+  void   SetNodeAndNext( const SMDS_MeshNode* n, PolyVertex& v, size_t index );
+  void   GetTriaNodes( const SMDS_MeshNode** nodes, size_t* nodeIndices) const;
+  double TriaArea() const;
+  bool   IsInsideTria( const PolyVertex* v );
+  PolyVertex* Delete();
+
+  struct Compare // compare PolyVertex'es by node
+  {
+    bool operator()(const PolyVertex* a, const PolyVertex* b) const
+    {
+      return ( a->_nxyz.Node() <  b->_nxyz.Node() );
+    }
+  };
+  // set of PolyVertex sorted by mesh node
+  typedef boost::container::flat_set< PolyVertex*, Compare > PVSet;
+};
+
+struct Triangulate::Data
+{
+  std::vector< PolyVertex > _pv;
+  std::vector< size_t >     _nodeIndex;
+  PolyVertex::PVSet         _uniqueNodePV;
+};
 
 struct Triangulate::Optimizer
 {
-  std::vector< TNodeSet > _nodeUsage; // inclusions of a node in triangles
+  std::vector< TriaNodeSet > _nodeUsage; // inclusions of a node in triangles
 
   //================================================================================
   /*!
@@ -107,19 +141,19 @@ struct Triangulate::Optimizer
         size_t i2 = iTria + ( i + 1 ) % 3;
         size_t ind1 = nodeIndices[ i1 ]; // node index in points
         size_t ind2 = nodeIndices[ i2 ];
-        TNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
+        TriaNodeSet & usage1 = _nodeUsage[ ind1 ]; // triangles using a node
-        TNodeSet & usage2 = _nodeUsage[ ind2 ];
+        TriaNodeSet & usage2 = _nodeUsage[ ind2 ];
         if ( usage1.size() < 2 ||
              usage2.size() < 2 )
           continue;
 
         // look for another triangle using two nodes
-        TNodeSet::iterator usIt1 = usage1.begin();
+        TriaNodeSet::iterator usIt1 = usage1.begin();
         for ( ; usIt1 != usage1.end(); ++usIt1 )
         {
           if ( usIt1->_triaIndex == iTria )
             continue; // current triangle
-          TNodeSet::iterator usIt2 = usage2.find( *usIt1 );
+          TriaNodeSet::iterator usIt2 = usage2.find( *usIt1 );
           if ( usIt2 == usage2.end() )
             continue; // no common _triaIndex in two usages
 
@@ -138,13 +172,13 @@ struct Triangulate::Optimizer
           // swap edge by modifying nodeIndices
 
           nodeIndices[ i2 ] = ind4;
-          _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
           _nodeUsage[ ind4 ].insert({ iTria, i2 - iTria });
+          _nodeUsage[ ind2 ].erase ({ iTria, i2 - iTria });
 
           i1 = usIt1->Index();
           nodeIndices[ i1 ] = ind3;
-          _nodeUsage[ ind1 ].erase ( *usIt1 );
           _nodeUsage[ ind3 ].insert( *usIt1 );
+          _nodeUsage[ ind1 ].erase ( *usIt1 );
 
           --i; // to re-check a current edge
           badness1 = badness3;
@@ -170,16 +204,16 @@ struct Triangulate::Optimizer
                          std::vector< PolyVertex > & points,
                          bool                        checkArea = false )
   {
-    //if ( checkArea )
+    if ( checkArea )
     {
       points[ i2 ]._prev = & points[ i1 ];
       points[ i2 ]._next = & points[ i3 ];
       double a = points[ i2 ].TriaArea();
-      if ( a < 0 )
+      // if ( a < 0 )
-        return std::numeric_limits<double>::max();
+      //   return std::numeric_limits<double>::max();
-      return 1. / a;
+      // return 1. / a;
 
-      if ( points[ i2 ].TriaArea() < 0 )
+      if ( a < 0 )
         return 2;
     }
     const gp_XY & p1 = points[ i1 ]._xy;
@@ -311,12 +345,28 @@ bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v )
 bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
                                const size_t                        nbNodes)
 {
+  std::vector< PolyVertex >&    _pv = _data->_pv;
+  std::vector< size_t >& _nodeIndex = _data->_nodeIndex;
+  PolyVertex::PVSet&  _uniqueNodePV = _data->_uniqueNodePV;
+
   // connect nodes into a ring
   _pv.resize( nbNodes );
   for ( size_t i = 1; i < nbNodes; ++i )
     _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i], i-1 );
   _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0], nbNodes-1 );
 
+  // assure correctness of PolyVertex::_index as a node can encounter more than once
+  // within a polygon boundary
+  if ( _optimizer && nbNodes > 4 )
+  {
+    _uniqueNodePV.clear();
+    for ( size_t i = 0; i < nbNodes; ++i )
+    {
+      PolyVertex::PVSet::iterator pv = _uniqueNodePV.insert( &_pv[i] ).first;
+      _pv[i]._index = (*pv)->_index;
+    }
+  }
+
   // get a polygon normal
   gp_XYZ normal(0,0,0), p0,v01,v02;
   p0  = _pv[0]._nxyz;
@@ -342,11 +392,16 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
     _pv[i]._xy.SetY( axes.YDirection().XYZ() * p );
   }
 
+  // compute minimal triangle area
+  double sumArea = 0;
+  for ( size_t i = 0; i < nbNodes; ++i )
+    sumArea += _pv[i].TriaArea();
+  const double minArea = 1e-6 * sumArea / ( nbNodes - 2 );
+
   // in a loop, find triangles with positive area and having no vertices inside
   int iN = 0, nbTria = nbNodes - 2;
   nodes.resize( nbTria * 3 );
   _nodeIndex.resize( nbTria * 3 );
-  const double minArea = 1e-6;
   PolyVertex* v = &_pv[0], *vi;
   int nbVertices = nbNodes, nbBadTria = 0, isGoodTria;
   while ( nbBadTria < nbVertices )
@@ -430,6 +485,7 @@ bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes,
 
 Triangulate::Triangulate( bool optimize ): _optimizer(0)
 {
+  _data = new Data;
   if ( optimize )
     _optimizer = new Optimizer;
 }
@@ -442,6 +498,7 @@ Triangulate::Triangulate( bool optimize ): _optimizer(0)
 
 Triangulate::~Triangulate()
 {
+  delete _data;
   delete _optimizer;
   _optimizer = 0;
 }

src/SMESH_I/SMESH_Gen_i.cxx

@@ -1993,14 +1993,16 @@ CORBA::Boolean SMESH_Gen_i::Compute( SMESH::SMESH_Mesh_ptr theMesh,
 void SMESH_Gen_i::CancelCompute( SMESH::SMESH_Mesh_ptr theMesh,
                                  GEOM::GEOM_Object_ptr theShapeObject )
 {
-  SMESH_Mesh_i* meshServant = dynamic_cast<SMESH_Mesh_i*>( GetServant( theMesh ).in() );
+  if ( SMESH_Mesh_i* meshServant = dynamic_cast<SMESH_Mesh_i*>( GetServant( theMesh ).in() ))
-  ::SMESH_Mesh& myLocMesh = meshServant->GetImpl();
+  {
-  TopoDS_Shape myLocShape;
+    ::SMESH_Mesh& myLocMesh = meshServant->GetImpl();
-  if(theMesh->HasShapeToMesh())
+    TopoDS_Shape myLocShape;
-    myLocShape = GeomObjectToShape( theShapeObject );
+    if(theMesh->HasShapeToMesh())
-  else
+      myLocShape = GeomObjectToShape( theShapeObject );
-    myLocShape = SMESH_Mesh::PseudoShape();
+    else
-  myGen.CancelCompute( myLocMesh, myLocShape);
+      myLocShape = SMESH_Mesh::PseudoShape();
+    myGen.CancelCompute( myLocMesh, myLocShape);
+  }
 }
 
 //=============================================================================

src/SMESH_SWIG/smeshBuilder.py

@@ -6300,7 +6300,7 @@ class Mesh(metaclass = MeshMeta):
                 SubMeshOrGroup = [ obj.GetMesh() ]
                 break
             if isinstance( obj, int ):
-                SubMeshOrGroup = self.GetIDSource( SubMeshOrGroup, SMESH.NODE )
+                SubMeshOrGroup = [ self.GetIDSource( SubMeshOrGroup, SMESH.NODE )]
                 unRegister.set( SubMeshOrGroup )
                 break
 

src/StdMeshers/StdMeshers_Projection_2D.cxx

@@ -1694,7 +1694,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
   if ( !projDone || is1DComputed )
     // ----------------------------------------------------------------
     // The mapper can create distorted faces by placing nodes out of the FACE
-    // boundary, also bad face can be created if EDGEs already discretized
+    // boundary, also bad faces can be created if EDGEs already discretized
     // --> fix bad faces by smoothing
     // ----------------------------------------------------------------
     if ( helper.IsDistorted2D( tgtSubMesh, /*checkUV=*/false, &helper ))