Salome HOME
54355: 'Compute' button is absent for 'Number of the double nodes' value in 'Mesh...
[modules/smesh.git] / src / SMESHUtils / SMESH_MAT2d.cxx
index 1c89a7e0e853be737b1e2cfeb574dac4f3e132ec..8e3eaeb6481d4b44f0c31cd2c705641b30f5664d 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -47,6 +47,7 @@
 #include <TopoDS_Wire.hxx>
 
 #ifdef _DEBUG_
+//#define _MYDEBUG_
 #include "SMESH_File.hxx"
 #include "SMESH_Comment.hxx"
 #endif
@@ -76,6 +77,8 @@ namespace
 
     size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
     bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
+    bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
+                                                           Abs( _b - v->y() ) < 1. ); }
   };
   // -------------------------------------------------------------------------------------
 
@@ -90,8 +93,11 @@ namespace
     list< const TVDEdge* > _edges; // MA edges in CCW order within _cell
 
     InSegment( InPoint * p0, InPoint * p1, size_t iE)
-      : _p0(p0), _p1(p1), _geomEdgeInd(iE) {}
-    InSegment() : _p0(0), _p1(0), _geomEdgeInd(0) {}
+      : _p0(p0), _p1(p1), _geomEdgeInd(iE), _cell(0) {}
+    InSegment() : _p0(0), _p1(0), _geomEdgeInd(0), _cell(0) {}
+
+    const InPoint& point0() const { return *_p0; }
+    const InPoint& point1() const { return *_p1; }
 
     inline bool isConnected( const TVDEdge* edge );
 
@@ -105,10 +111,12 @@ namespace
   // check  if a TVDEdge begins at my end or ends at my start
   inline bool InSegment::isConnected( const TVDEdge* edge )
   {
-    return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
-             Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
-            (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
-             Abs( edge->vertex1()->y() - _p0->_b ) < 1.  ));
+    return (( edge->vertex0() && edge->vertex1() )
+            &&
+            ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
+              Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+             (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+              Abs( edge->vertex1()->y() - _p0->_b ) < 1.  )));
   }
 
   // check if a MA TVDEdge is outside of a domain
@@ -147,14 +155,17 @@ namespace
   // }
 
   // -------------------------------------------------------------------------------------
-#ifdef _DEBUG_
+#ifdef _MYDEBUG_
   // writes segments into a txt file readable by voronoi_visualizer
   void inSegmentsToFile( vector< InSegment>& inSegments)
   {
     if ( inSegments.size() > 1000 )
       return;
     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
+    const char*     user = getenv("USER");
+    if ( !user || strcmp( user, "eap" )) return;
     SMESH_File file(fileName, false );
+    file.remove();
     file.openForWriting();
     SMESH_Comment text;
     text << "0\n"; // nb points
@@ -180,7 +191,7 @@ namespace
     if ( !edge->vertex1() )
       cout << ") -> ( INF, INF";
     else
-      cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+      cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
     cout << ")\t cell=" << edge->cell()
          << " iBnd=" << edge->color()
          << " twin=" << edge->twin()
@@ -209,9 +220,9 @@ namespace
     }
   }
 #else
-  void inSegmentsToFile( vector< InSegment>& inSegments) {}
-  void dumpEdge( const TVDEdge* edge ) {}
-  void dumpCell( const TVDCell* cell ) {}
+  #define inSegmentsToFile(arg) {}
+  //void dumpEdge( const TVDEdge* edge ) {}
+  //void dumpCell( const TVDCell* cell ) {}
 #endif
 }
 // -------------------------------------------------------------------------------------
@@ -252,7 +263,9 @@ namespace boost {
 
 namespace
 {
-  const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+  const int    theNoBrachID = 0;
+  double       theScale[2]; // scale used in bndSegsToMesh()
+  const size_t theNoEdgeID = std::numeric_limits<size_t>::max() / 1000;
 
   // -------------------------------------------------------------------------------------
   /*!
@@ -270,17 +283,18 @@ namespace
   };
   // -------------------------------------------------------------------------------------
   /*!
-   * \brief A segment on EDGE, used to create BndPoints
+   * \brief Segment of EDGE, used to create BndPoints
    */
   struct BndSeg
   {
     InSegment*       _inSeg;
     const TVDEdge*   _edge;
     double           _uLast;
+    BndSeg*          _prev; // previous BndSeg in FACE boundary
     int              _branchID; // negative ID means reverse direction
 
     BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
-      _inSeg(seg), _edge(edge), _uLast(u), _branchID( theNoBrachID ) {}
+      _inSeg(seg), _edge(edge), _uLast(u), _prev(0), _branchID( theNoBrachID ) {}
 
     void setIndexToEdge( size_t id )
     {
@@ -291,29 +305,50 @@ namespace
 
     size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
 
-    void setBranch( int branchID, vector< BndSeg >& bndSegs )
+    static BndSeg* getBndSegOfEdge( const TVDEdge*              edge,
+                                    vector< vector< BndSeg > >& bndSegsPerEdge )
     {
-      _branchID = branchID;
-
-      if ( _edge ) // pass branch to an opposite BndSeg
+      BndSeg* seg = 0;
+      if ( edge )
       {
-        size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
-        if ( oppSegIndex < bndSegs.size() && bndSegs[ oppSegIndex ]._branchID == theNoBrachID )
-          bndSegs[ oppSegIndex ]._branchID = -branchID;
+        size_t oppSegIndex  = SMESH_MAT2d::Branch::getBndSegment( edge );
+        size_t oppEdgeIndex = SMESH_MAT2d::Branch::getGeomEdge  ( edge );
+        if ( oppEdgeIndex < bndSegsPerEdge.size() &&
+             oppSegIndex < bndSegsPerEdge[ oppEdgeIndex ].size() )
+        {
+          seg = & bndSegsPerEdge[ oppEdgeIndex ][ oppSegIndex ];
+        }
       }
+      return seg;
     }
-    bool hasOppositeEdge( const size_t noEdgeID )
+
+    void setBranch( int branchID, vector< vector< BndSeg > >& bndSegsPerEdge )
+    {
+      _branchID = branchID;
+
+      // pass branch to an opposite BndSeg
+      if ( _edge )
+        if ( BndSeg* oppSeg = getBndSegOfEdge( _edge->twin(), bndSegsPerEdge ))
+        {
+          if ( oppSeg->_branchID == theNoBrachID )
+            oppSeg->_branchID = -branchID;
+        }
+    }
+    bool hasOppositeEdge()
     {
       if ( !_edge ) return false;
-      return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != noEdgeID );
+      return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != theNoEdgeID );
     }
 
-    // check a next segment in CW order
+    // check a next segment in CCW order
     bool isSameBranch( const BndSeg& seg2 )
     {
       if ( !_edge || !seg2._edge )
         return true;
 
+      if ( _edge->twin() == seg2._edge )
+        return true;
+
       const TVDCell* cell1 = this->_edge->twin()->cell();
       const TVDCell* cell2 = seg2. _edge->twin()->cell();
       if ( cell1 == cell2 )
@@ -337,19 +372,110 @@ namespace
       else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
       {
         if ( edgeMedium1->twin() == edgeMedium2 &&
-             SMESH_MAT2d::Branch::getBndSegment( edgeMedium1 ) ==
-             SMESH_MAT2d::Branch::getBndSegment( edgeMedium2 ))
+             SMESH_MAT2d::Branch::getGeomEdge( edgeMedium1 ) ==
+             SMESH_MAT2d::Branch::getGeomEdge( edgeMedium2 ))
           // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
           return true;
       }
 
       return false;
     }
+  }; // struct BndSeg
+
+  // -------------------------------------------------------------------------------------
+  /*!
+   * \brief Iterator 
+   */
+  struct BranchIterator
+  {
+    int                                 _i, _size;
+    const std::vector<const TVDEdge*> & _edges;
+    bool                                _closed;
+
+    BranchIterator(const std::vector<const TVDEdge*> & edges, int i )
+      :_i( i ), _size( edges.size() ), _edges( edges )
+    {
+      _closed = ( edges[0]->vertex1() == edges.back()->vertex0() || // closed branch
+                  edges[0]->vertex0() == edges.back()->vertex1() );
+    }
+    const TVDEdge* operator++() { ++_i; return edge(); }
+    const TVDEdge* operator--() { --_i; return edge(); }
+    bool operator<( const BranchIterator& other ) { return _i < other._i; }
+    BranchIterator& operator=( const BranchIterator& other ) { _i = other._i; return *this; }
+    void set(int i) { _i = i; }
+    int  index() const { return _i; }
+    int  indexMod() const { return ( _i + _size ) % _size; }
+    const TVDEdge* edge() const {
+      return _closed ? _edges[ indexMod() ] : ( _i < 0 || _i >= _size ) ? 0 : _edges[ _i ];
+    }
+    const TVDEdge* edgePrev() { --_i; const TVDEdge* e = edge(); ++_i; return e; }
+    const TVDEdge* edgeNext() { ++_i; const TVDEdge* e = edge(); --_i; return e; }
   };
 
   //================================================================================
   /*!
-   * \brief Computes length of of TVDEdge
+   * \brief debug: to visually check found MA edges
+   */
+  //================================================================================
+
+  void bndSegsToMesh( const vector< vector< BndSeg > >& bndSegsPerEdge )
+  {
+#ifdef _MYDEBUG_
+    if ( !getenv("bndSegsToMesh")) return;
+    map< const TVDVertex *, int > v2Node;
+    map< const TVDVertex *, int >::iterator v2n;
+    set< const TVDEdge* > addedEdges;
+
+    const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
+    SMESH_File file(fileName, false );
+    file.remove();
+    file.openForWriting();
+    SMESH_Comment text;
+    text << "import salome, SMESH\n";
+    text << "salome.salome_init()\n";
+    text << "from salome.smesh import smeshBuilder\n";
+    text << "smesh = smeshBuilder.New()\n";
+    text << "m=smesh.Mesh()\n";
+    for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
+    {
+      const vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
+      {
+        if ( !bndSegs[i]._edge )
+          text << "# E=" << iE << " i=" << i << " NULL edge\n";
+        else if ( !bndSegs[i]._edge->vertex0() ||
+                  !bndSegs[i]._edge->vertex1() )
+          text << "# E=" << iE << " i=" << i << " INFINITE edge\n";
+        else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
+                  addedEdges.insert( bndSegs[i]._edge->twin() ).second )
+        {
+          v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
+          size_t n0 = v2n->second;
+          if ( n0 == v2Node.size() )
+            text << "n" << n0 << " = m.AddNode( "
+                 << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
+                 << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
+
+          v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
+          size_t n1 = v2n->second;
+          if ( n1 == v2Node.size() )
+            text << "n" << n1 << " = m.AddNode( "
+                 << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
+                 << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
+
+          text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
+        }
+      }
+    }
+    text << "\n";
+    file.write( text.c_str(), text.size() );
+    cout << fileName << endl;
+#endif
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes length of a TVDEdge
    */
   //================================================================================
 
@@ -522,10 +648,18 @@ namespace
     for ( size_t iE = 0; iE < edges.size(); ++iE )
     {
       size_t iE2 = (iE+1) % edges.size();
-      if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared ))
-        continue;
+      if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared )) // FACE with several WIREs?
+        for ( size_t i = 1; i < edges.size(); ++i )
+        {
+          iE2 = (iE2+1) % edges.size();
+          if ( iE != iE2 &&
+               TopExp::CommonVertex( edges[iE], edges[iE2], vShared ) &&
+               vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
+            break;
+        }
       if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
-        return false;
+        continue;
+        //return false;
       vector< UVU > & points1 = uvuVec[ iE ];
       vector< UVU > & points2 = uvuVec[ iE2 ];
       gp_Pnt2d & uv1 = points1.back() ._uv;
@@ -536,7 +670,7 @@ namespace
     // get scale to have the same 2d proportions as in 3d
     computeProportionScale( face, uvBox, scale );
 
-    // make scale to have coordinates precise enough when converted to int
+    // make 'scale' such that to have coordinates precise enough when converted to int
 
     gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
     uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
@@ -546,7 +680,7 @@ namespace
     double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )),
                        Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) };
     int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1;
-    const double precision = 1e-5;
+    const double precision = Min( 1e-5, minSegLen * 1e-2 );
     double preciScale = Min( vMax[iMax] / precision,
                              std::numeric_limits<int>::max() / vMax[iMax] );
     preciScale /= scale[iMax];
@@ -577,6 +711,8 @@ namespace
         {
           inPoints[ iP++ ] = points[i-1].getInPoint( scale );
           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+          if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
+            return false; // too short segment
         }
       }
     }
@@ -590,9 +726,15 @@ namespace
         {
           inPoints[ iP++ ] = points[i].getInPoint( scale );
           inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+          if ( inPoints[ iP-2 ] == inPoints[ iP-1 ])
+            return false; // too short segment
         }
       }
     }
+    // debug
+    theScale[0] = scale[0];
+    theScale[1] = scale[1];
+
     return true;
   }
 
@@ -602,22 +744,25 @@ namespace
    */
   //================================================================================
 
-  void updateJoinedBranch( vector< const TVDEdge* > & branchEdges,
-                           const size_t               newID,
-                           vector< BndSeg > &         bndSegs,
-                           const bool                 reverse)
+  void updateJoinedBranch( vector< const TVDEdge* > &   branchEdges,
+                           const size_t                 newID,
+                           vector< vector< BndSeg > > & bndSegs,
+                           const bool                   reverse)
   {
+    BndSeg *seg1, *seg2;
     if ( reverse )
     {
       for ( size_t i = 0; i < branchEdges.size(); ++i )
       {
-        size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
-        size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
-        bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
-        bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
-        bndSegs[ seg1 ]._branchID *= -newID;
-        bndSegs[ seg2 ]._branchID *= -newID;
-        branchEdges[i] = branchEdges[i]->twin();
+        if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i],         bndSegs )) &&
+            ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
+        {
+          seg1->_branchID /= seg1->branchID();
+          seg2->_branchID /= seg2->branchID();
+          seg1->_branchID *= -newID;
+          seg2->_branchID *= -newID;
+          branchEdges[i] = branchEdges[i]->twin();
+        }
       }
       std::reverse( branchEdges.begin(), branchEdges.end() );
     }
@@ -625,12 +770,14 @@ namespace
     {
       for ( size_t i = 0; i < branchEdges.size(); ++i )
       {
-        size_t seg1 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i] );
-        size_t seg2 = SMESH_MAT2d::Branch::getBndSegment( branchEdges[i]->twin() );
-        bndSegs[ seg1 ]._branchID /= bndSegs[ seg1 ].branchID();
-        bndSegs[ seg2 ]._branchID /= bndSegs[ seg2 ].branchID();
-        bndSegs[ seg1 ]._branchID *= newID;
-        bndSegs[ seg2 ]._branchID *= newID;
+        if (( seg1 = BndSeg::getBndSegOfEdge( branchEdges[i],         bndSegs )) &&
+            ( seg2 = BndSeg::getBndSegOfEdge( branchEdges[i]->twin(), bndSegs )))
+        {
+          seg1->_branchID /= seg1->branchID();
+          seg2->_branchID /= seg2->branchID();
+          seg1->_branchID *= newID;
+          seg2->_branchID *= newID;
+        }
       }
     }
   }
@@ -655,12 +802,16 @@ namespace
                vector< const SMESH_MAT2d::BranchEnd* >& branchPnt,
                SMESH_MAT2d::Boundary&                   boundary )
   {
-    const size_t noEdgeID = inSegments.size() + 1; // ID of non-existent geom EDGE
-
-    // Associate MA cells with inSegments
+    // Associate MA cells with geom EDGEs
     for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
     {
       const TVDCell* cell = &(*it);
+      if ( cell->is_degenerate() )
+      {
+        std::cerr << "SMESH_MAT2d: encounter degenerate voronoi_cell. Invalid input data?"
+                  << std::endl;
+        return;
+      }
       if ( cell->contains_segment() )
       {
         InSegment& seg = inSegments[ cell->source_index() ];
@@ -669,7 +820,7 @@ namespace
       }
       else
       {
-        InSegment::setGeomEdgeToCell( cell, noEdgeID );
+        InSegment::setGeomEdgeToCell( cell, theNoEdgeID );
       }
     }
 
@@ -729,13 +880,16 @@ namespace
           continue;
         inPntChecked[ pInd ] = true;
 
-        const TVDEdge* edge =  // a TVDEdge passing through an end of inSeg
-          is2nd ? maEdges.front()->prev() : maEdges.back()->next();
-        while ( true )
+        const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
+        if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
+          continue;
+        const TVDEdge* edge =  // a secondary TVDEdge connecting inPnt and maE
+          is2nd ? maE->prev() : maE->next();
+        while ( inSeg.isConnected( edge ))
         {
           if ( edge->is_primary() ) break; // this should not happen
           const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
-          if ( inSeg.getGeomEdge( edge2->cell() ) != noEdgeID )
+          if ( inSeg.getGeomEdge( edge2->cell() ) != theNoEdgeID )
             break; // cell of an InSegment
           bool hasInfinite = false;
           list< const TVDEdge* > pointEdges;
@@ -769,61 +923,96 @@ namespace
     if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ )
     {
       inPntChecked[0] = false; // do not use the 1st point twice
-      //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), noEdgeID );
+      //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), theNoEdgeID );
       inPoints[0]._edges.clear();
     }
 
     // Divide InSegment's into BndSeg's (so that each BndSeg corresponds to one MA edge)
 
-    vector< BndSeg > bndSegs;
-    bndSegs.reserve( inSegments.size() * 3 );
-
-    list< const TVDEdge* >::reverse_iterator e;
-    for ( size_t i = 0; i < inSegments.size(); ++i )
+    vector< vector< BndSeg > > bndSegsPerEdge( boundary.nbEdges() ); // all BndSeg's
     {
-      InSegment& inSeg = inSegments[i];
+      vector< BndSeg > bndSegs; // bndSeg's of a current EDGE
+      size_t prevGeomEdge = theNoEdgeID;
 
-      // segments around 1st concave point
-      size_t ip0 = inSeg._p0->index( inPoints );
-      if ( inPntChecked[ ip0 ] )
-        for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
-          bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
-      inPntChecked[ ip0 ] = false;
-
-      // segments of InSegment's
-      const size_t nbMaEdges = inSeg._edges.size();
-      switch ( nbMaEdges ) {
-      case 0: // "around" circle center
-        bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
-      case 1:
-        bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
-      default:
-        gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
-                        inSeg._p1->_b - inSeg._p0->_b );
-        const double inSegLen2 = inSegDir.SquareModulus();
-        e = inSeg._edges.rbegin();
-        for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
+      list< const TVDEdge* >::reverse_iterator e;
+      for ( size_t i = 0; i < inSegments.size(); ++i )
+      {
+        InSegment& inSeg = inSegments[i];
+
+        if ( inSeg._geomEdgeInd != prevGeomEdge )
         {
-          gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
-                      (*e)->vertex0()->y() - inSeg._p0->_b );
-          double r = toMA * inSegDir / inSegLen2;
-          double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
-          bndSegs.push_back( BndSeg( &inSeg, *e, u ));
+          if ( !bndSegs.empty() )
+            bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
+          prevGeomEdge = inSeg._geomEdgeInd;
         }
-        bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
-      }
-      // segments around 2nd concave point
-      size_t ip1 = inSeg._p1->index( inPoints );
-      if ( inPntChecked[ ip1 ] )
-        for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
+
+        // segments around 1st concave point
+        size_t ip0 = inSeg._p0->index( inPoints );
+        if ( inPntChecked[ ip0 ] )
+          for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e )
+            bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param ));
+        inPntChecked[ ip0 ] = false;
+
+        // segments of InSegment's
+        const size_t nbMaEdges = inSeg._edges.size();
+        switch ( nbMaEdges ) {
+        case 0: // "around" circle center
+          bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break;
+        case 1:
+          bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break;
+        default:
+          gp_XY inSegDir( inSeg._p1->_a - inSeg._p0->_a,
+                          inSeg._p1->_b - inSeg._p0->_b );
+          const double inSegLen2 = inSegDir.SquareModulus();
+          e = inSeg._edges.rbegin();
+          for ( size_t iE = 1; iE < nbMaEdges; ++e, ++iE )
+          {
+            gp_XY toMA( (*e)->vertex0()->x() - inSeg._p0->_a,
+                        (*e)->vertex0()->y() - inSeg._p0->_b );
+            double r = toMA * inSegDir / inSegLen2;
+            double u = r * inSeg._p1->_param + ( 1. - r ) * inSeg._p0->_param;
+            bndSegs.push_back( BndSeg( &inSeg, *e, u ));
+          }
           bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
-      inPntChecked[ ip1 ] = false;
+        }
+        // segments around 2nd concave point
+        size_t ip1 = inSeg._p1->index( inPoints );
+        if ( inPntChecked[ ip1 ] )
+          for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e )
+            bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param ));
+        inPntChecked[ ip1 ] = false;
+      }
+      if ( !bndSegs.empty() )
+        bndSegsPerEdge[ prevGeomEdge ].swap( bndSegs );
     }
 
-    // make TVDEdge's know it's BndSeg to enable passing branchID to
-    // an opposite BndSeg in BndSeg::setBranch()
-    for ( size_t i = 0; i < bndSegs.size(); ++i )
-      bndSegs[i].setIndexToEdge( i );
+    // prepare to MA branch search
+    for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
+    {
+      // 1) make TVDEdge's know it's BndSeg to enable passing branchID to
+      // an opposite BndSeg in BndSeg::setBranch(); geom EDGE ID is known from TVDCell
+      // 2) connect bndSegs via BndSeg::_prev
+
+      vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
+      if ( bndSegs.empty() ) continue;
+
+      for ( size_t i = 1; i < bndSegs.size(); ++i )
+      {
+        bndSegs[i]._prev = & bndSegs[i-1];
+        bndSegs[i].setIndexToEdge( i );
+      }
+      // look for the last bndSeg of previous EDGE to set bndSegs[0]._prev
+      const InPoint& p0 = bndSegs[0]._inSeg->point0();
+      for ( size_t iE2 = 0; iE2 < bndSegsPerEdge.size(); ++iE2 )
+        if ( p0 == bndSegsPerEdge[ iE2 ].back()._inSeg->point1() )
+        {
+          bndSegs[0]._prev = & bndSegsPerEdge[ iE2 ].back();
+          break;
+        }
+      bndSegs[0].setIndexToEdge( 0 );
+    }
+
+    bndSegsToMesh( bndSegsPerEdge ); // debug: visually check found MA edges
 
 
     // Find TVDEdge's of Branches and associate them with bndSegs
@@ -834,76 +1023,117 @@ namespace
     map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
 
     int branchID = 1; // we code orientation as branchID sign
-    branchEdges.resize( branchID + 1 );
+    branchEdges.resize( branchID );
 
-    size_t i1st = 0;
-    while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
-      ++i1st;
-    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and the opposite bndSeg
-    branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
+    vector< std::pair< int, const TVDVertex* > > branchesToCheckEnd;
 
-    for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
+    for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
     {
-      if ( bndSegs[i].branchID() )
+      vector< BndSeg >& bndSegs = bndSegsPerEdge[ iE ];
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
       {
-        branchID = bndSegs[i]._branchID; // with sign
+        if ( bndSegs[i].branchID() )
+        {
+          if ( bndSegs[i]._prev &&
+               bndSegs[i]._branchID == -bndSegs[i]._prev->_branchID &&
+               bndSegs[i]._edge )
+          {
+            SMESH_MAT2d::BranchEndType type =
+              ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
+                SMESH_MAT2d::BE_ON_VERTEX :
+                SMESH_MAT2d::BE_END );
+            endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
+          }
+          continue;
+        }
+        if ( !bndSegs[i]._prev &&
+             !bndSegs[i].hasOppositeEdge() )
+          continue;
 
-        if ( bndSegs[i]._branchID == -bndSegs[i-1]._branchID &&
-             bndSegs[i]._edge )
+        if ( !bndSegs[i]._prev ||
+             !bndSegs[i]._prev->isSameBranch( bndSegs[i] ))
         {
-          SMESH_MAT2d::BranchEndType type =
-            ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ?
-              SMESH_MAT2d::BE_ON_VERTEX :
-              SMESH_MAT2d::BE_END );
-          endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type ));
+          branchEdges.resize(( branchID = branchEdges.size()) + 1 );
+          if ( bndSegs[i]._edge && bndSegs[i]._prev )
+          {
+            endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_BRANCH_POINT ));
+            if ( bndSegs[i]._prev->_branchID < 0 )
+              // 0023404: a branch-point is inside a branch
+              branchesToCheckEnd.push_back( make_pair( bndSegs[i]._prev->branchID(),
+                                                       bndSegs[i]._edge->vertex1() ));
+          }
         }
-        continue;
-      }
-      if ( !bndSegs[i-1].isSameBranch( bndSegs[i] ))
-      {
-        branchEdges.resize(( branchID = branchEdges.size()) + 1 );
-        if ( bndSegs[i]._edge )
-          endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
-                                     SMESH_MAT2d::BE_BRANCH_POINT ));
+        else if ( bndSegs[i]._prev->_branchID )
+        {
+          branchID = bndSegs[i]._prev->_branchID;  // with sign
+        }
+        else if ( bndSegs[i]._edge ) // 1st bndSeg of a WIRE
+        {
+          branchEdges.resize(( branchID = branchEdges.size()) + 1 );
+          if ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ))
+          {
+            if ( bndSegs[i]._inSeg->point0() == bndSegs[i]._edge->vertex1() )
+              endType.insert( make_pair( bndSegs[i]._edge->vertex1(), SMESH_MAT2d::BE_ON_VERTEX ));
+            else
+              endType.insert( make_pair( bndSegs[i]._edge->vertex0(), SMESH_MAT2d::BE_ON_VERTEX ));
+          }
+        }
+
+        bndSegs[i].setBranch( branchID, bndSegsPerEdge ); // set to i-th and to the opposite bndSeg
+        if ( bndSegs[i].hasOppositeEdge() )
+          branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
       }
-      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
-      if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
-        branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
     }
-    // define BranchEndType of the first TVDVertex
-    if ( bndSegs.front()._branchID == -bndSegs.back()._branchID )
+
+    if ( !ignoreCorners && !branchesToCheckEnd.empty() )
     {
-      if ( bndSegs[0]._edge )
+      // split branches having branch-point inside
+      // (a branch-point was not detected since another branch is joined at the opposite side)
+      for ( size_t i = 0; i < branchesToCheckEnd.size(); ++i )
       {
-        SMESH_MAT2d::BranchEndType type =
-          ( bndSegs[0]._inSeg->isConnected( bndSegs[0]._edge ) ?
-            SMESH_MAT2d::BE_ON_VERTEX :
-            SMESH_MAT2d::BE_END );
-        endType.insert( make_pair( bndSegs[0]._edge->vertex1(), type ));
-      }
-      else if ( bndSegs.back()._edge )
-      {
-        SMESH_MAT2d::BranchEndType type =
-          ( bndSegs.back()._inSeg->isConnected( bndSegs.back()._edge ) ?
-            SMESH_MAT2d::BE_ON_VERTEX :
-            SMESH_MAT2d::BE_END );
-        endType.insert( make_pair( bndSegs.back()._edge->vertex0(), type ));
+        vector<const TVDEdge*> & branch = branchEdges[ branchesToCheckEnd[i].first ];
+        const TVDVertex*    branchPoint = branchesToCheckEnd[i].second;
+        if ( branch.front()->vertex1() == branchPoint ||
+             branch.back ()->vertex0() == branchPoint )
+          continue; // OK - branchPoint is at a branch end
+
+        // find a MA edge where another branch begins
+        size_t iE;
+        for ( iE = 0; iE < branch.size(); ++iE )
+          if ( branch[iE]->vertex1() == branchPoint )
+            break;
+        if ( iE < branch.size() )
+        {
+          // split the branch
+          branchEdges.resize(( branchID = branchEdges.size()) + 1 );
+          vector<const TVDEdge*> & branch2 = branchEdges[ branchID ];
+          branch2.assign( branch.begin()+iE, branch.end() );
+          branch.resize( iE );
+          for ( iE = 0; iE < branch2.size(); ++iE )
+            if ( BndSeg* bs = BndSeg::getBndSegOfEdge( branch2[iE], bndSegsPerEdge ))
+              bs->setBranch( branchID, bndSegsPerEdge );
+        }
       }
     }
+
     // join the 1st and the last branch edges if it is the same branch
-    if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
-         bndSegs.back().isSameBranch( bndSegs.front() ))
-    {
-      vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
-      vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID()  ];
-      br1.insert( br1.begin(), br2.begin(), br2.end() );
-      br2.clear();
-    }
+    // if ( bndSegs.back().branchID() != bndSegs.front().branchID() &&
+    //      bndSegs.back().isSameBranch( bndSegs.front() ))
+    // {
+    //   vector<const TVDEdge*> & br1 = branchEdges[ bndSegs.front().branchID() ];
+    //   vector<const TVDEdge*> & br2 = branchEdges[ bndSegs.back().branchID()  ];
+    //   br1.insert( br1.begin(), br2.begin(), br2.end() );
+    //   br2.clear();
+    // }
 
-    // remove branches ending at BE_ON_VERTEX
+    // remove branches ending at BE_ON_VERTEX and BE_END
 
     vector<bool> isBranchRemoved( branchEdges.size(), false );
 
+    std::set< SMESH_MAT2d::BranchEndType > endTypeToRm;
+    endTypeToRm.insert( SMESH_MAT2d::BE_ON_VERTEX );
+    endTypeToRm.insert( SMESH_MAT2d::BE_END );
+
     if ( ignoreCorners && branchEdges.size() > 2 && !branchEdges[2].empty() )
     {
       // find branches to remove
@@ -915,10 +1145,10 @@ namespace
         const TVDVertex* v0 = branchEdges[iB][0]->vertex1();
         const TVDVertex* v1 = branchEdges[iB].back()->vertex0();
         v2et = endType.find( v0 );
-        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+        if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
           isBranchRemoved[ iB ] = true;
         v2et = endType.find( v1 );
-        if ( v2et != endType.end() && v2et->second == SMESH_MAT2d::BE_ON_VERTEX )
+        if ( v2et != endType.end() && endTypeToRm.count( v2et->second ))
           isBranchRemoved[ iB ] = true;
       }
       // try to join not removed branches into one
@@ -937,34 +1167,38 @@ namespace
         if ( !v0 && !v1 )
           continue;
 
-        size_t iBrToJoin = 0;
-        for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
+        for ( int isV0 = 0; isV0 < 2; ++isV0 )
         {
-          if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
-            continue;
-          const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
-          const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
-          if ( v0 == v02 || v0 == v12 || v1 == v02 || v1 == v12 )
+          const TVDVertex* v = isV0 ? v0 : v1;
+          size_t iBrToJoin = 0;
+          for ( size_t iB2 = 1; iB2 < branchEdges.size(); ++iB2 )
           {
-            if ( iBrToJoin > 0 )
+            if ( branchEdges[iB2].empty() || isBranchRemoved[iB2] || iB == iB2 )
+              continue;
+            const TVDVertex* v02 = branchEdges[iB2][0]->vertex1();
+            const TVDVertex* v12 = branchEdges[iB2].back()->vertex0();
+            if ( v == v02 || v == v12 )
             {
-              iBrToJoin = 0;
-              break; // more than 2 not removed branches meat at a TVDVertex
+              if ( iBrToJoin > 0 )
+              {
+                iBrToJoin = 0;
+                break; // more than 2 not removed branches meat at a TVDVertex
+              }
+              iBrToJoin = iB2;
             }
-            iBrToJoin = iB2;
           }
-        }
-        if ( iBrToJoin > 0 )
-        {
-          vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
-          const TVDVertex* v02 = branch[0]->vertex1();
-          const TVDVertex* v12 = branch.back()->vertex0();
-          updateJoinedBranch( branch, iB, bndSegs, /*reverse=*/(v0 == v02 || v1 == v12 ));
-          if ( v0 == v02 || v0 == v12 )
-            branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
-          else
-            branchEdges[iB].insert( branchEdges[iB].end(),   branch.begin(), branch.end() );
-          branch.clear();
+          if ( iBrToJoin > 0 )
+          {
+            vector<const TVDEdge*>& branch = branchEdges[ iBrToJoin ];
+            const TVDVertex* v02 = branch[0]->vertex1();
+            const TVDVertex* v12 = branch.back()->vertex0();
+            updateJoinedBranch( branch, iB, bndSegsPerEdge, /*reverse=*/(v0 == v02 || v1 == v12 ));
+            if ( v0 == v02 || v0 == v12 )
+              branchEdges[iB].insert( branchEdges[iB].begin(), branch.begin(), branch.end() );
+            else
+              branchEdges[iB].insert( branchEdges[iB].end(),   branch.begin(), branch.end() );
+            branch.clear();
+          }
         }
       } // loop on branchEdges
     } // if ( ignoreCorners )
@@ -991,36 +1225,31 @@ namespace
 
     // Fill in BndPoints of each EDGE of the boundary
 
-    size_t iSeg = 0;
+    //size_t iSeg = 0;
     int edgeInd = -1, dInd = 0;
-    while ( iSeg < bndSegs.size() )
+    for ( size_t iE = 0; iE < bndSegsPerEdge.size(); ++iE )
     {
-      const size_t                geomID = bndSegs[ iSeg ].geomEdge();
-      SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( geomID );
-
-      size_t nbSegs = 0;
-      for ( size_t i = iSeg; i < bndSegs.size() && geomID == bndSegs[ i ].geomEdge(); ++i )
-        ++nbSegs;
-      size_t iSegEnd = iSeg + nbSegs;
+      vector< BndSeg >&          bndSegs = bndSegsPerEdge[ iE ];
+      SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( iE );
 
       // make TVDEdge know an index of bndSegs within BndPoints
-      for ( size_t i = iSeg; i < iSegEnd; ++i )
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
         if ( bndSegs[i]._edge )
-          SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge );
+          SMESH_MAT2d::Branch::setBndSegment( i, bndSegs[i]._edge );
 
       // parameters on EDGE
 
-      bndPoints._params.reserve( nbSegs + 1 );
-      bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param );
+      bndPoints._params.reserve( bndSegs.size() + 1 );
+      bndPoints._params.push_back( bndSegs[ 0 ]._inSeg->_p0->_param );
 
-      for ( size_t i = iSeg; i < iSegEnd; ++i )
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
         bndPoints._params.push_back( bndSegs[ i ]._uLast );
 
       // MA edges
 
-      bndPoints._maEdges.reserve( nbSegs );
+      bndPoints._maEdges.reserve( bndSegs.size() );
 
-      for ( size_t i = iSeg; i < iSegEnd; ++i )
+      for ( size_t i = 0; i < bndSegs.size(); ++i )
       {
         const size_t              brID = bndSegs[ i ].branchID();
         const SMESH_MAT2d::Branch*  br = branchByID[ brID ];
@@ -1044,7 +1273,7 @@ namespace
             else // bndSegs[ i ]._branchID > 0
             {
               dInd = +1;
-              for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
+              for ( edgeInd = 0; edgeInd < (int)branchEdges[ brID ].size(); ++edgeInd )
                 if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
                   break;
             }
@@ -1063,10 +1292,7 @@ namespace
         bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
 
       } // loop on bndSegs of an EDGE
-
-      iSeg = iSegEnd;
-
-    } // loop on all bndSegs
+    } // loop on all bndSegs to construct Boundary
 
     // Initialize branches
 
@@ -1218,13 +1444,22 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
     while ( points._params[i+1] < u ) ++i;
   }
 
+  if ( points._params[i] == points._params[i+1] ) // coincident points at some end
+  {
+    int di = ( points._params[0] == points._params[i] ) ? +1 : -1;
+    while ( points._params[i] == points._params[i+1] )
+      i += di;
+    if ( i < 0 || i+1 >= (int)points._params.size() )
+      i = 0;
+  }
+
   double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
 
   if ( !points._maEdges[ i ].second ) // no branch at the EDGE end, look for a closest branch
   {
-    if ( i < points._maEdges.size() / 2 ) // near 1st point
+    if ( i < (int)points._maEdges.size() / 2 ) // near 1st point
     {
-      while ( i < points._maEdges.size()-1 && !points._maEdges[ i ].second )
+      while ( i < (int)points._maEdges.size()-1 && !points._maEdges[ i ].second )
         ++i;
       edgeParam = edgeReverse;
     }
@@ -1245,6 +1480,21 @@ bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Returns a BranchPoint corresponding to a given BoundaryPoint on a geom EDGE
+ *  \param [in] bp - the BoundaryPoint
+ *  \param [out] p - the found BranchPoint
+ *  \return bool - is OK
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getBranchPoint( const BoundaryPoint& bp,
+                                            BranchPoint&         p ) const
+{
+  return getBranchPoint( bp._edgeIndex, bp._param, p );
+}
+
 //================================================================================
 /*!
  * \brief Check if a given boundary segment is a null-length segment on a concave
@@ -1279,7 +1529,7 @@ bool SMESH_MAT2d::Boundary::moveToClosestEdgeEnd( BoundaryPoint& bp ) const
     return false;
 
   const BndPoints& points = _pointsPerEdge[ bp._edgeIndex ];
-  if ( bp._param - points._params[0] < points._params.back() - bp._param )
+  if ( Abs( bp._param - points._params[0]) < Abs( points._params.back() - bp._param ))
     bp._param = points._params[0];
   else 
     bp._param = points._params.back();
@@ -1331,9 +1581,9 @@ Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) cons
  */
 //================================================================================
 
-void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                maEdges,
-                                const Boundary*                        boundary,
-                                map< const TVDVertex*, BranchEndType > endType )
+void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                 maEdges,
+                                const Boundary*                         boundary,
+                                map< const TVDVertex*, BranchEndType >& endType )
 {
   if ( maEdges.empty() ) return;
 
@@ -1512,7 +1762,7 @@ bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
   if ( p._iEdge > _params.size()-1 )
     return false;
   if ( p._iEdge == _params.size()-1 )
-    return u = 1.;
+    return ( u = 1. );
 
   u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
         _params[ p._iEdge+1 ] * p._edgeParam );
@@ -1585,7 +1835,7 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
                                                    std::vector< BranchPoint >&   divPoints,
                                                    const vector<const TVDEdge*>& maEdges,
                                                    const vector<const TVDEdge*>& maEdgesTwin,
-                                                   size_t &                    i) const
+                                                   int &                         i) const
 {
   // if there is a concave vertex between EDGEs
   // then position of a dividing BranchPoint is undefined, it is somewhere
@@ -1597,11 +1847,13 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
   BranchPoint divisionPnt;
   divisionPnt._branch = this;
 
+  BranchIterator iCur( maEdges, i );
+
   size_t ie1 = getGeomEdge( maEdges    [i] );
   size_t ie2 = getGeomEdge( maEdgesTwin[i] );
 
-  size_t iSeg1  = getBndSegment( maEdges[ i-1 ] );
-  size_t iSeg2  = getBndSegment( maEdges[ i ] );
+  size_t iSeg1  = getBndSegment( iCur.edgePrev() );
+  size_t iSeg2  = getBndSegment( iCur.edge() );
   bool isConcaPrev = _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 );
   bool isConcaNext = _boundary->isConcaveSegment( ie1,             iSeg2 );
   if ( !isConcaNext && !isConcaPrev )
@@ -1609,46 +1861,48 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
 
   bool isConcaveV = false;
 
-  int iPrev = i-1, iNext = i;
+  const TVDEdge* maE;
+  BranchIterator iPrev( maEdges, i ), iNext( maEdges, i );
+   --iPrev;
   if ( isConcaNext ) // all null-length segments follow
   {
     // look for a VERTEX of the opposite EDGE
-    ++iNext; // end of null-length segments
-    while ( iNext < maEdges.size() )
+    // iNext - next after all null-length segments
+    while (( maE = ++iNext ))
     {
-      iSeg2 = getBndSegment( maEdges[ iNext ] );
-      if ( _boundary->isConcaveSegment( ie1, iSeg2 ))
-        ++iNext;
-      else
+      iSeg2 = getBndSegment( maE );
+      if ( !_boundary->isConcaveSegment( ie1, iSeg2 ))
         break;
     }
     bool vertexFound = false;
-    for ( size_t iE = i+1; iE < iNext; ++iE )
+    for ( ++iCur; iCur < iNext; ++iCur )
     {
-      ie2 = getGeomEdge( maEdgesTwin[iE] );
+      ie2 = getGeomEdge( maEdgesTwin[ iCur.indexMod() ] );
       if ( ie2 != edgeIDs2.back() )
       {
         // opposite VERTEX found
-        divisionPnt._iEdge = iE;
+        divisionPnt._iEdge = iCur.indexMod();
         divisionPnt._edgeParam = 0;
         divPoints.push_back( divisionPnt );
         edgeIDs1.push_back( ie1 );
         edgeIDs2.push_back( ie2 );
         vertexFound = true;
-        }
+      }
     }
     if ( vertexFound )
     {
-      iPrev = i = --iNext; // not to add a BP in the moddle
+      --iNext;
+      iPrev = iNext; // not to add a BP in the moddle
+      i = iNext.indexMod();
       isConcaveV = true;
     }
   }
   else if ( isConcaPrev )
   {
     // all null-length segments passed, find their beginning
-    while ( iPrev-1 >= 0 )
+    while (( maE = iPrev.edgePrev() ))
     {
-      iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
+      iSeg1 = getBndSegment( maE );
       if ( _boundary->isConcaveSegment( edgeIDs1.back(), iSeg1 ))
         --iPrev;
       else
@@ -1656,17 +1910,18 @@ bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&
     }
   }
 
-  if ( iPrev < i-1 || iNext > i )
+  if ( iPrev.index() < i-1 || iNext.index() > i )
   {
     // no VERTEX on the opposite EDGE, put the Branch Point in the middle
-    double par1 = _params[ iPrev+1 ], par2 = _params[ iNext ];
+    divisionPnt._iEdge = iPrev.indexMod();
+    ++iPrev;
+    double   par1 = _params[ iPrev.indexMod() ], par2 = _params[ iNext.indexMod() ];
     double midPar = 0.5 * ( par1 + par2 );
-    divisionPnt._iEdge = iPrev;
-    while ( _params[ divisionPnt._iEdge + 1 ] < midPar )
-      ++divisionPnt._iEdge;
+    for ( ; _params[ iPrev.indexMod() ] < midPar; ++iPrev )
+      divisionPnt._iEdge = iPrev.indexMod();
     divisionPnt._edgeParam =
-      ( _params[ divisionPnt._iEdge + 1 ] - midPar ) /
-      ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] );
+      ( _params[ iPrev.indexMod() ] - midPar ) /
+      ( _params[ iPrev.indexMod() ] - _params[ divisionPnt._iEdge ] );
     divPoints.push_back( divisionPnt );
     isConcaveV = true;
   }
@@ -1694,50 +1949,74 @@ void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edge
   edgeIDs2.clear();
   divPoints.clear();
 
-  edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
-  edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
-
   std::vector<const TVDEdge*> twins( _maEdges.size() );
   for ( size_t i = 0; i < _maEdges.size(); ++i )
     twins[i] = _maEdges[i]->twin();
 
+  BranchIterator maIter ( _maEdges, 0 );
+  BranchIterator twIter ( twins, 0 );
   // size_t lastConcaE1 = _boundary.nbEdges();
   // size_t lastConcaE2 = _boundary.nbEdges();
 
+  // if ( maIter._closed ) // closed branch
+  // {
+  //   edgeIDs1.push_back( getGeomEdge( _maEdges.back() ));
+  //   edgeIDs2.push_back( getGeomEdge( _maEdges.back()->twin() ));
+  // }
+  // else
+  {
+    edgeIDs1.push_back( getGeomEdge( maIter.edge() ));
+    edgeIDs2.push_back( getGeomEdge( twIter.edge() ));
+  }
+
   BranchPoint divisionPnt;
   divisionPnt._branch = this;
 
-  for ( size_t i = 0; i < _maEdges.size(); ++i )
+  for ( ++maIter, ++twIter; maIter.index() < (int)_maEdges.size(); ++maIter, ++twIter )
   {
-    size_t ie1 = getGeomEdge( _maEdges[i] );
-    size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
-    
-    if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
+    size_t ie1 = getGeomEdge( maIter.edge() );
+    size_t ie2 = getGeomEdge( twIter.edge() );
+
+    bool otherE1 = ( edgeIDs1.back() != ie1 );
+    bool otherE2 = ( edgeIDs2.back() != ie2 );
+
+    if ( !otherE1 && !otherE2 && maIter._closed )
+    {
+      int iSegPrev1 = getBndSegment( maIter.edgePrev() );
+      int iSegCur1  = getBndSegment( maIter.edge() );
+      otherE1 = Abs( iSegPrev1 - iSegCur1 ) != 1;
+      int iSegPrev2 = getBndSegment( twIter.edgePrev() );
+      int iSegCur2  = getBndSegment( twIter.edge() );
+      otherE2 = Abs( iSegPrev2 - iSegCur2 ) != 1;
+    }
+
+    if ( otherE1 || otherE2 )
     {
       bool isConcaveV = false;
-      if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 )
+      if ( otherE1 && !otherE2 )
       {
-        isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i );
+        isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints,
+                                              _maEdges, twins, maIter._i );
       }
-      if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 )
+      if ( !otherE1 && otherE2 )
       {
-        isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i );
+        isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints,
+                                              twins, _maEdges, maIter._i );
       }
 
       if ( isConcaveV )
       {
-        ie1 = getGeomEdge( _maEdges[i] );
-        ie2 = getGeomEdge( _maEdges[i]->twin() );
+        ie1 = getGeomEdge( maIter.edge() );
+        ie2 = getGeomEdge( twIter.edge() );
       }
-      if (( !isConcaveV ) ||
-          ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ))
+      if ( !isConcaveV || otherE1 || otherE2 )
       {
         edgeIDs1.push_back( ie1 );
         edgeIDs2.push_back( ie2 );
       }
       if ( divPoints.size() < edgeIDs1.size() - 1 )
       {
-        divisionPnt._iEdge = i;
+        divisionPnt._iEdge = maIter.index();
         divisionPnt._edgeParam = 0;
         divPoints.push_back( divisionPnt );
       }