diff --git a/doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png b/doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png
new file mode 100644
index 000000000..b02ceabcd
Binary files /dev/null and b/doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png differ
diff --git a/doc/salome/gui/SMESH/images/quad_from_ma_mesh.png b/doc/salome/gui/SMESH/images/quad_from_ma_mesh.png
new file mode 100644
index 000000000..f233cc640
Binary files /dev/null and b/doc/salome/gui/SMESH/images/quad_from_ma_mesh.png differ
diff --git a/doc/salome/gui/SMESH/input/basic_meshing_algos.doc b/doc/salome/gui/SMESH/input/basic_meshing_algos.doc
index f494926d1..b82b67f89 100644
--- a/doc/salome/gui/SMESH/input/basic_meshing_algos.doc
+++ b/doc/salome/gui/SMESH/input/basic_meshing_algos.doc
@@ -58,7 +58,8 @@ without geometrical objects.
 
 There is also a number of more specific algorithms:
 <ul>
-<li>\subpage prism_3d_algo_page "for meshing prismatic shapes"</li>
+<li>\subpage prism_3d_algo_page "for meshing prismatic 3D shapes"</li>
+<li>\subpage quad_from_ma_algo_page "for meshing faces with sinuous borders"</li>
 <li>\subpage projection_algos_page "for meshing by projection of another mesh"</li>
 <li>\subpage import_algos_page "for meshing by importing elements from another mesh"</li>
 <li>\subpage radial_prism_algo_page "for meshing geometrical objects with cavities"</li>
diff --git a/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc b/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc
new file mode 100644
index 000000000..976783d40
--- /dev/null
+++ b/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc
@@ -0,0 +1,30 @@
+/*!
+
+\page quad_from_ma_algo_page Medial Axis Projection Quadrangle meshing algorithm
+
+Medial Axis Projection algorithm can be used for meshing faces with
+sinuous borders and having channel-like shape, for which is it
+difficult to define 1D hypotheses so that generated quadrangles to be
+of good shape.
+
+\image html quad_from_ma_mesh.png "A mesh of a river model"
+
+The algorithm assures good shape of quadrangles by constructing Medial
+Axis between sinuous borders of the face and using it to
+discretize the borders.
+
+\image html quad_from_ma_medial_axis.png "Media Axis between two blue sinuous borders"
+
+The Medial Axis is used in two ways:
+<ol>
+<li>If there is a sub-mesh on either sinuous border, then the nodes of
+  this border are mapped to the opposite border via the Medial
+  Axis.</li>
+<li> If there is no sub-meshes on the sinuous borders, then a part of
+  the Medial Axis that can be mapped to both borders is discretized
+  using a hypothesis assigned to the face or its ancestor shapes,
+  and the division points are mapped from the Medial Axis to the both
+  borders.</li>
+</ol>
+
+*/
diff --git a/idl/SMESH_BasicHypothesis.idl b/idl/SMESH_BasicHypothesis.idl
index 274365bf2..60f603caa 100644
--- a/idl/SMESH_BasicHypothesis.idl
+++ b/idl/SMESH_BasicHypothesis.idl
@@ -1090,6 +1090,13 @@ module StdMeshers
   {
   };
 
+  /*!
+   * StdMeshers_QuadFromMedialAxis_1D2D: interface of "Quadrangle (Medial Axis Projection)" algorithm
+   */
+  interface StdMeshers_QuadFromMedialAxis_1D2D : SMESH::SMESH_2D_Algo
+  {
+  };
+
   /*!
    * StdMeshers_Hexa_3D: interface of "Hexahedron (i,j,k)" algorithm
    */
diff --git a/resources/StdMeshers.xml.in b/resources/StdMeshers.xml.in
index 350b9ba8c..4150b42dc 100644
--- a/resources/StdMeshers.xml.in
+++ b/resources/StdMeshers.xml.in
@@ -311,6 +311,19 @@
       </python-wrap>
     </algorithm>
 
+    <algorithm type     ="QuadFromMedialAxis_1D2D"
+               label-id ="Quadrangle (Medial Axis Projection)"
+               icon-id  ="mesh_algo_quad.png"
+               opt-hypos="ViscousLayers2D"
+               input    ="EDGE"
+               output   ="QUAD"
+               dim      ="2">
+      <python-wrap>
+        <algo>QuadFromMedialAxis_1D2D=Quadrangle(algo=smeshBuilder.QUAD_MA_PROJ)</algo>
+        <hypo>ViscousLayers2D=ViscousLayers2D(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetIgnoreEdges())</hypo>
+      </python-wrap>
+    </algorithm>
+
     <algorithm type     ="Hexa_3D"
                label-id ="Hexahedron (i,j,k)"
                icon-id  ="mesh_algo_hexa.png"
diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx
index ab4c5832c..7537bf8a8 100644
--- a/src/SMESH/SMESH_MesherHelper.cxx
+++ b/src/SMESH/SMESH_MesherHelper.cxx
@@ -32,6 +32,7 @@
 #include "SMDS_IteratorOnIterators.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESH_Block.hxx"
+#include "SMESH_HypoFilter.hxx"
 #include "SMESH_MeshAlgos.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
@@ -3206,6 +3207,28 @@ TopAbs_ShapeEnum SMESH_MesherHelper::GetGroupType(const TopoDS_Shape& group,
   return TopAbs_SHAPE;
 }
 
+//================================================================================
+/*!
+ * \brief Returns a shape, to which a hypothesis used to mesh a given shape is assigned
+ *  \param [in] hyp - the hypothesis
+ *  \param [in] shape - the shape, for meshing which the \a hyp is used
+ *  \param [in] mesh - the mesh
+ *  \return TopoDS_Shape - the shape the \a hyp is assigned to
+ */
+//================================================================================
+
+TopoDS_Shape SMESH_MesherHelper::GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp,
+                                                       const TopoDS_Shape&        shape,
+                                                       SMESH_Mesh*                mesh)
+{
+  const SMESH_Hypothesis* h = static_cast<const SMESH_Hypothesis*>( hyp );
+  SMESH_HypoFilter hypFilter( SMESH_HypoFilter::Is( h ));
+
+  TopoDS_Shape shapeOfHyp;
+  mesh->GetHypothesis( shape, hypFilter, /*checkAncestors=*/true, &shapeOfHyp );
+  return shapeOfHyp;
+}
+
 //=======================================================================
 //function : IsQuadraticMesh
 //purpose  : Check mesh without geometry for: if all elements on this shape are quadratic,
diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx
index 9b17d6bb9..05173ae06 100644
--- a/src/SMESH/SMESH_MesherHelper.hxx
+++ b/src/SMESH/SMESH_MesherHelper.hxx
@@ -236,6 +236,10 @@ class SMESH_EXPORT SMESH_MesherHelper
   static TopAbs_ShapeEnum GetGroupType(const TopoDS_Shape& group,
                                        const bool          avoidCompound=false);
 
+  static TopoDS_Shape GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp,
+                                            const TopoDS_Shape&        shape,
+                                            SMESH_Mesh*                mesh);
+
 
 public:
   // ---------- PUBLIC INSTANCE METHODS ----------
@@ -243,7 +247,9 @@ public:
   // constructor
   SMESH_MesherHelper(SMESH_Mesh& theMesh);
 
-  SMESH_Mesh* GetMesh() const { return myMesh; }
+  SMESH_Gen*    GetGen() const { return GetMesh()->GetGen(); }
+    
+  SMESH_Mesh*   GetMesh() const { return myMesh; }
     
   SMESHDS_Mesh* GetMeshDS() const { return GetMesh()->GetMeshDS(); }
     
diff --git a/src/SMESHUtils/CMakeLists.txt b/src/SMESHUtils/CMakeLists.txt
index 486454e74..8c2f3217d 100644
--- a/src/SMESHUtils/CMakeLists.txt
+++ b/src/SMESHUtils/CMakeLists.txt
@@ -63,6 +63,7 @@ SET(SMESHUtils_HEADERS
   SMESH_Utils.hxx
   SMESH_TryCatch.hxx
   SMESH_MeshAlgos.hxx
+  SMESH_MAT2d.hxx
 )
 
 # --- sources ---
@@ -76,6 +77,7 @@ SET(SMESHUtils_SOURCES
   SMESH_TryCatch.cxx
   SMESH_File.cxx
   SMESH_MeshAlgos.cxx
+  SMESH_MAT2d.cxx
 )
 
 # --- rules ---
diff --git a/src/SMESHUtils/SMESH_MAT2d.cxx b/src/SMESHUtils/SMESH_MAT2d.cxx
new file mode 100644
index 000000000..390d337b8
--- /dev/null
+++ b/src/SMESHUtils/SMESH_MAT2d.cxx
@@ -0,0 +1,1571 @@
+// Copyright (C) 2007-2015  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
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : SMESH_MAT2d.cxx
+// Created   : Thu May 28 17:49:53 2015
+// Author    : Edward AGAPOV (eap)
+
+#include "SMESH_MAT2d.hxx"
+
+#include <list>
+
+#include <BRepAdaptor_CompCurve.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <Bnd_B2d.hxx>
+//#include <GCPnts_AbscissaPoint.hxx>
+#include <GCPnts_TangentialDeflection.hxx>
+// #include <GCPnts_UniformAbscissa.hxx>
+// #include <GCPnts_UniformDeflection.hxx>
+#include <Geom2d_Curve.hxx>
+//#include <GeomAdaptor_Curve.hxx>
+#include <Geom2dAdaptor_Curve.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Surface.hxx>
+#include <TopExp.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <TopoDS_Wire.hxx>
+
+#ifdef _DEBUG_
+#include "SMESH_File.hxx"
+#include "SMESH_Comment.hxx"
+#endif
+
+using namespace std;
+using boost::polygon::x;
+using boost::polygon::y;
+using SMESH_MAT2d::TVD;
+using SMESH_MAT2d::TVDEdge;
+using SMESH_MAT2d::TVDCell;
+using SMESH_MAT2d::TVDVertex;
+
+namespace
+{
+  // Input data for construct_voronoi()
+  // -------------------------------------------------------------------------------------
+
+  struct InPoint
+  {
+    int _a, _b;
+    double _param;
+    InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {}
+    InPoint() : _a(0), _b(0), _param(0) {}
+
+    // working data
+    list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order
+
+    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; }
+  };
+  // -------------------------------------------------------------------------------------
+
+  struct InSegment
+  {
+    InPoint * _p0;
+    InPoint * _p1;
+
+    // working data
+    size_t                 _geomEdgeInd; // EDGE index within the FACE
+    const TVDCell*         _cell;
+    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) {}
+
+    inline bool isConnected( const TVDEdge* edge );
+
+    inline bool isExternal( const TVDEdge* edge );
+
+    static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); }
+
+    static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); }
+  };
+
+  // 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.  ));
+  }
+
+  // check if a MA TVDEdge is outside of a domain
+  inline bool InSegment::isExternal( const TVDEdge* edge )
+  {
+    double dot = // x1*x2 + y1*y2;   (x1,y1) - internal normal of InSegment
+      ( _p0->_b - _p1->_b ) * ( 0.5 * ( edge->vertex0()->x() + edge->vertex1()->x() ) - _p0->_a ) +
+      ( _p1->_a - _p0->_a ) * ( 0.5 * ( edge->vertex0()->y() + edge->vertex1()->y() ) - _p0->_b );
+    return dot < 0.;
+  }
+
+  // // -------------------------------------------------------------------------------------
+  // const size_t theExternMA = 111; // to mark external MA edges
+
+  // bool isExternal( const TVDEdge* edge )
+  // {
+  //   return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA );
+  // }
+
+  // // mark external MA edges
+  // void markExternalEdges( const TVDEdge* edge )
+  // {
+  //   if ( isExternal( edge ))
+  //     return;
+  //   SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge );
+  //   SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() );
+  //   if ( edge->is_primary() && edge->vertex1() )
+  //   {
+  //     const TVDVertex * v = edge->vertex1();
+  //     edge = v->incident_edge();
+  //     do {
+  //       markExternalEdges( edge );
+  //       edge = edge->rot_next();
+  //     } while ( edge != v->incident_edge() );
+  //   }
+  // }
+
+  // -------------------------------------------------------------------------------------
+#ifdef _DEBUG_
+  // 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";
+    SMESH_File file(fileName, false );
+    file.openForWriting();
+    SMESH_Comment text;
+    text << "0\n"; // nb points
+    text << inSegments.size() << "\n"; // nb segments
+    for ( size_t i = 0; i < inSegments.size(); ++i )
+    {
+      text << inSegments[i]._p0->_a << " "
+           << inSegments[i]._p0->_b << " "
+           << inSegments[i]._p1->_a << " "
+           << inSegments[i]._p1->_b << "\n";
+    }
+    text << "\n";
+    file.write( text.c_str(), text.size() );
+    cout << "Write " << fileName << endl;
+  }
+  void dumpEdge( const TVDEdge* edge )
+  {
+    cout << "*Edge_" << edge;
+    if ( !edge->vertex0() )
+      cout << " ( INF, INF";
+    else
+      cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y();
+    if ( !edge->vertex1() )
+      cout << ") -> ( INF, INF";
+    else
+      cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+    cout << ")\t cell=" << edge->cell()
+         << " iBnd=" << edge->color()
+         << " twin=" << edge->twin()
+         << " twin_cell=" << edge->twin()->cell()
+         << " prev=" << edge->prev() << " next=" << edge->next()
+         << ( edge->is_primary() ? " MA " : " SCND" )
+         << ( edge->is_linear() ? " LIN " : " CURV" )
+         << endl;
+  }
+  void dumpCell( const TVDCell* cell )
+  {
+    cout << "**Cell_" << cell << " GEOM=" << cell->color() << " ";
+    cout << ( cell->contains_segment() ? " SEG " : " PNT " );
+    if ( cell-> is_degenerate() )
+      cout << " degen ";
+    else
+    {
+      cout << endl;
+      const TVDEdge* edge = cell->incident_edge();
+      size_t i = 0;
+      do {
+        edge = edge->next();
+        cout << "   - " << ++i << " ";
+        dumpEdge( edge );
+      } while (edge != cell->incident_edge());
+    }
+  }
+#else
+  void inSegmentsToFile( vector< InSegment>& inSegments) {}
+  void dumpEdge( const TVDedge* edge ) {}
+  void dumpCell( const TVDCell* cell ) {}
+#endif
+}
+// -------------------------------------------------------------------------------------
+
+namespace boost {
+  namespace polygon {
+
+    template <>
+    struct geometry_concept<InPoint> {
+      typedef point_concept type;
+    };
+    template <>
+    struct point_traits<InPoint> {
+      typedef int coordinate_type;
+
+      static inline coordinate_type get(const InPoint& point, orientation_2d orient) {
+        return (orient == HORIZONTAL) ? point._a : point._b;
+      }
+    };
+
+    template <>
+    struct geometry_concept<InSegment> {
+      typedef segment_concept type;
+    };
+
+    template <>
+    struct segment_traits<InSegment> {
+      typedef int coordinate_type;
+      typedef InPoint point_type;
+
+      static inline point_type get(const InSegment& segment, direction_1d dir) {
+        return *(dir.to_int() ? segment._p1 : segment._p0);
+      }
+    };
+  }  // namespace polygon
+} // namespace boost
+  // -------------------------------------------------------------------------------------
+
+namespace
+{
+  const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+
+  // -------------------------------------------------------------------------------------
+  /*!
+   * \brief Intermediate DS to create InPoint's
+   */
+  struct UVU
+  {
+    gp_Pnt2d _uv;
+    double   _u;
+    UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {}
+    InPoint getInPoint( double scale[2] )
+    {
+      return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u );
+    }
+  };
+  // -------------------------------------------------------------------------------------
+  /*!
+   * \brief A segment on EDGE, used to create BndPoints
+   */
+  struct BndSeg
+  {
+    InSegment*       _inSeg;
+    const TVDEdge*   _edge;
+    double           _uLast;
+    int              _branchID; // negative ID means reverse direction
+
+    BndSeg( InSegment* seg, const TVDEdge* edge, double u ):
+      _inSeg(seg), _edge(edge), _uLast(u), _branchID( theNoBrachID ) {}
+
+    void setIndexToEdge( size_t id )
+    {
+      SMESH_MAT2d::Branch::setBndSegment( id, _edge );
+    }
+
+    int branchID() const { return Abs( _branchID ); }
+
+    size_t geomEdge() const { return _inSeg->_geomEdgeInd; }
+
+    void setBranch( int branchID, vector< BndSeg >& bndSegs )
+    {
+      _branchID = branchID;
+
+      if ( _edge ) // pass branch to an opposite BndSeg
+      {
+        size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() );
+        if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ )
+          bndSegs[ oppSegIndex ]._branchID = -branchID;
+      }
+    }
+    bool hasOppositeEdge( const size_t noEdgeID )
+    {
+      if ( !_edge ) return false;
+      return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != noEdgeID );
+    }
+
+    // check a next segment in CW order
+    bool isSameBranch( const BndSeg& seg2 )
+    {
+      if ( !_edge || !seg2._edge )
+        return true;
+
+      const TVDCell* cell1 = this->_edge->twin()->cell();
+      const TVDCell* cell2 = seg2. _edge->twin()->cell();
+      if ( cell1 == cell2 )
+        return true;
+
+      const TVDEdge* edgeMedium1 = this->_edge->twin()->next();
+      const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev();
+
+      if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() )
+      {
+        if ( edgeMedium1->twin() == edgeMedium2 )
+          return true;
+        // edgeMedium's are edges whose twin()->cell is built on an end point of inSegment
+        // and is located between cell1 and cell2
+        if ( edgeMedium1->twin() == edgeMedium2->twin() ) // is this possible???
+          return true;
+        if ( edgeMedium1->twin() == edgeMedium2->twin()->next() &&
+             edgeMedium1->twin()->cell()->contains_point() )
+          return true;
+      }
+      else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() )
+      {
+        if ( edgeMedium1->twin() == edgeMedium2 &&
+             SMESH_MAT2d::Branch::getBndSegment( edgeMedium1 ) ==
+             SMESH_MAT2d::Branch::getBndSegment( edgeMedium2 ))
+          // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner
+          return true;
+      }
+
+      return false;
+    }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Computes length of of TVDEdge
+   */
+  //================================================================================
+
+  double length( const TVDEdge* edge )
+  {
+    gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(),
+             edge->vertex0()->y() - edge->vertex1()->y() );
+    return d.Modulus();
+  }
+
+  //================================================================================
+  /*!
+   * \brief Compute scale to have the same 2d proportions as in 3d
+   */
+  //================================================================================
+
+  void computeProportionScale( const TopoDS_Face& face,
+                               const Bnd_B2d&     uvBox,
+                               double             scale[2])
+  {
+    scale[0] = scale[1] = 1.;
+    if ( uvBox.IsVoid() ) return;
+
+    TopLoc_Location loc;
+    Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
+
+    const int nbDiv = 30;
+    gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
+    gp_XY uvMid = 0.5 * ( uvMin + uvMax );
+    double du = ( uvMax.X() - uvMin.X() ) / nbDiv;
+    double dv = ( uvMax.Y() - uvMin.Y() ) / nbDiv;
+
+    double uLen3d = 0, vLen3d = 0;
+    gp_Pnt uPrevP = surface->Value( uvMin.X(), uvMid.Y() );
+    gp_Pnt vPrevP = surface->Value( uvMid.X(), uvMin.Y() );
+    for (int i = 1; i <= nbDiv; i++)
+    {
+      double u = uvMin.X() + du * i;
+      double v = uvMin.Y() + dv * i;
+      gp_Pnt uP = surface->Value( u, uvMid.Y() );
+      gp_Pnt vP = surface->Value( uvMid.X(), v );
+      uLen3d += uP.Distance( uPrevP );
+      vLen3d += vP.Distance( vPrevP );
+      uPrevP = uP;
+      vPrevP = vP;
+    }
+    scale[0] = uLen3d / ( uvMax.X() - uvMin.X() );
+    scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Fill input data for construct_voronoi()
+   */
+  //================================================================================
+
+  bool makeInputData(const TopoDS_Face&                face,
+                     const std::vector< TopoDS_Edge >& edges,
+                     const double                      minSegLen,
+                     vector< InPoint >&                inPoints,
+                     vector< InSegment>&               inSegments,
+                     double                            scale[2])
+  {
+    const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization
+    TopLoc_Location loc;
+
+    // discretize the EDGEs to get 2d points and segments
+
+    vector< vector< UVU > > uvuVec( edges.size() );
+    Bnd_B2d uvBox;
+    for ( size_t iE = 0; iE < edges.size(); ++iE )
+    {
+      vector< UVU > & points = uvuVec[ iE ];
+
+      double f,l;
+      Handle(Geom_Curve)   c3d = BRep_Tool::Curve         ( edges[ iE ], loc, f, l );
+      Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface( edges[ iE ], face, f, l );
+      if ( c2d.IsNull() ) return false;
+
+      points.push_back( UVU( c2d->Value( f ), f ));
+      uvBox.Add( points.back()._uv );
+
+      Geom2dAdaptor_Curve c2dAdaptor (c2d, f,l );
+      double curDeflect = 0.3; //0.01;  //Curvature deflection
+      double angDeflect = 0.2; // 0.09; //Angular deflection
+
+      GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect);
+      // if ( discret.NbPoints() > 2 )
+      // {
+      //   cout << endl;
+      //   do
+      //   {
+      //     discret.Initialize( c2dAdaptor, 100, curDeflect );
+      //     cout << "C " << curDeflect << " " << discret.NbPoints() << endl;
+      //     curDeflect *= 1.5;
+      //   }
+      //   while ( discret.NbPoints() > 5 );
+      //   cout << endl;
+      //   do
+      //   {
+      //     discret.Initialize( c2dAdaptor, angDeflect, 100 );
+      //     cout << "A " << angDeflect << " " << discret.NbPoints() << endl;
+      //     angDeflect *= 1.5;
+      //   }
+      //   while ( discret.NbPoints() > 5 );
+      // }
+      gp_Pnt p, pPrev;
+      if ( !c3d.IsNull() )
+        pPrev = c3d->Value( f );
+      for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point
+      {
+        double u = discret.Parameter(i);
+        if ( !c3d.IsNull() )
+        {
+          p = c3d->Value( u );
+          int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef );
+          double dU = ( u - points.back()._u ) / nbDiv;
+          for ( int iD = 1; iD < nbDiv; ++iD )
+          {
+            double uD = points.back()._u + dU;
+            points.push_back( UVU( c2d->Value( uD ), uD ));
+          }
+          pPrev = p;
+        }
+        points.push_back( UVU( c2d->Value( u ), u ));
+        uvBox.Add( points.back()._uv );
+      }
+      // if ( !c3d.IsNull() )
+      // {
+      //   vector<double> params;
+      //   GeomAdaptor_Curve c3dAdaptor( c3d,f,l );
+      //   if ( useDefl )
+      //   {
+      //     const double deflection = minSegLen * 0.1;
+      //     GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true );
+      //     if ( !discret.IsDone() )
+      //       return false;
+      //     int nbP = discret.NbPoints();
+      //     for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
+      //       params.push_back( discret.Parameter(i) );
+      //   }
+      //   else
+      //   {
+      //     double   eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor );
+      //     int     nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef ));
+      //     double segLen = eLen / nbSeg;
+      //     GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l );
+      //     int nbP = Min( discret.NbPoints(), nbSeg + 1 );
+      //     for ( int i = 2; i < nbP; i++ ) // skip 1st and last points
+      //       params.push_back( discret.Parameter(i) );
+      //   }
+      //   for ( size_t i = 0; i < params.size(); ++i )
+      //   {
+      //     points.push_back( UVU( c2d->Value( params[i] ), params[i] ));
+      //     uvBox.Add( points.back()._uv );
+      //   }
+      // }
+      if ( points.size() < 2 )
+      {
+        points.push_back( UVU( c2d->Value( l ), l ));
+        uvBox.Add( points.back()._uv );
+      }
+      if ( edges[ iE ].Orientation() == TopAbs_REVERSED )
+        std::reverse( points.begin(), points.end() );
+    }
+
+    // make connected EDGEs have same UV at shared VERTEX
+    TopoDS_Vertex vShared;
+    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 ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true )))
+        return false;
+      vector< UVU > & points1 = uvuVec[ iE ];
+      vector< UVU > & points2 = uvuVec[ iE2 ];
+      gp_Pnt2d & uv1 = points1.back() ._uv;
+      gp_Pnt2d & uv2 = points2.front()._uv;
+      uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() );
+    }
+
+    // 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
+
+    gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax();
+    uvMin.ChangeCoord(1) = uvMin.X() * scale[0];
+    uvMin.ChangeCoord(2) = uvMin.Y() * scale[1];
+    uvMax.ChangeCoord(1) = uvMax.X() * scale[0];
+    uvMax.ChangeCoord(2) = uvMax.Y() * scale[1];
+    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;
+    double preciScale = Min( vMax[iMax] / precision,
+                             std::numeric_limits<int>::max() / vMax[iMax] );
+    preciScale /= scale[iMax];
+    double roundedScale = 10; // to ease debug
+    while ( roundedScale * 10 < preciScale )
+      roundedScale *= 10.;
+    scale[0] *= roundedScale;
+    scale[1] *= roundedScale;
+
+    // create input points and segments
+
+    inPoints.clear();
+    inSegments.clear();
+    size_t nbPnt = 0;
+    for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
+      nbPnt += uvuVec[ iE ].size();
+    inPoints.resize( nbPnt );
+    inSegments.reserve( nbPnt );
+
+    size_t iP = 0;
+    if ( face.Orientation() == TopAbs_REVERSED )
+    {
+      for ( int iE = uvuVec.size()-1; iE >= 0; --iE )
+      {
+        vector< UVU > & points = uvuVec[ iE ];
+        inPoints[ iP++ ] = points.back().getInPoint( scale );
+        for ( size_t i = points.size()-1; i >= 1; --i )
+        {
+          inPoints[ iP++ ] = points[i-1].getInPoint( scale );
+          inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+        }
+      }
+    }
+    else
+    {
+      for ( size_t iE = 0; iE < uvuVec.size(); ++iE )
+      {
+        vector< UVU > & points = uvuVec[ iE ];
+        inPoints[ iP++ ] = points[0].getInPoint( scale );
+        for ( size_t i = 1; i < points.size(); ++i )
+        {
+          inPoints[ iP++ ] = points[i].getInPoint( scale );
+          inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE ));
+        }
+      }
+    }
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Create MA branches and FACE boundary data
+   *  \param [in] vd - voronoi diagram of \a inSegments
+   *  \param [in] inPoints - FACE boundary points
+   *  \param [in,out] inSegments - FACE boundary segments
+   *  \param [out] branch - MA branches to fill
+   *  \param [out] branchEnd - ends of MA branches to fill
+   *  \param [out] boundary - FACE boundary to fill
+   */
+  //================================================================================
+
+  void makeMA( const TVD&                               vd,
+               vector< InPoint >&                       inPoints,
+               vector< InSegment > &                    inSegments,
+               vector< SMESH_MAT2d::Branch >&           branch,
+               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
+    for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it)
+    {
+      const TVDCell* cell = &(*it);
+      if ( cell->contains_segment() )
+      {
+        InSegment& seg = inSegments[ cell->source_index() ];
+        seg._cell = cell;
+        seg.setGeomEdgeToCell( cell, seg._geomEdgeInd );
+      }
+      else
+      {
+        InSegment::setGeomEdgeToCell( cell, noEdgeID );
+      }
+    }
+
+    vector< bool > inPntChecked( inPoints.size(), false );
+
+    // Find MA edges of each inSegment
+
+    for ( size_t i = 0; i < inSegments.size(); ++i )
+    {
+      InSegment& inSeg = inSegments[i];
+
+      // get edges around the cell lying on MA
+      bool hasSecondary = false;
+      const TVDEdge* edge = inSeg._cell->incident_edge();
+      do {
+        edge = edge->next(); // Returns the CCW next edge within the cell.
+        if ( edge->is_primary() && !inSeg.isExternal( edge ) )
+          inSeg._edges.push_back( edge ); // edge equidistant from two InSegments
+        else
+          hasSecondary = true;
+      } while (edge != inSeg._cell->incident_edge());
+
+      // there can be several continuous MA edges but maEdges can begin in the middle of
+      // a chain of continuous MA edges. Make the chain continuous.
+      list< const TVDEdge* >& maEdges = inSeg._edges;
+      if ( maEdges.empty() )
+        continue;
+      if ( hasSecondary )
+        while ( maEdges.back()->next() == maEdges.front() )
+          maEdges.splice( maEdges.end(), maEdges, maEdges.begin() );
+
+      // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE
+      list< const TVDEdge* >::iterator e = maEdges.begin();
+      while ( e != maEdges.end() )
+      {
+        const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge
+        size_t         geoE2 = InSegment::getGeomEdge( cell2 );
+        bool        toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e ));
+        if ( toRemove )
+          e = maEdges.erase( e );
+        else
+          ++e;
+      }
+      if ( maEdges.empty() )
+        continue;
+
+      // add MA edges corresponding to concave InPoints
+      for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg
+      {
+        InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 );
+        size_t    pInd = inPnt.index( inPoints );
+        if ( inPntChecked[ pInd ] )
+          continue;
+        if ( pInd > 0 &&
+             inPntChecked[ pInd-1 ] &&
+             inPoints[ pInd-1 ] == inPnt )
+          continue;
+        inPntChecked[ pInd ] = true;
+
+        const TVDEdge* edge =  // a TVDEdge passing through an end of inSeg
+          is2nd ? maEdges.front()->prev() : maEdges.back()->next();
+        while ( true )
+        {
+          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 )
+            break; // cell of an InSegment
+          bool hasInfinite = false;
+          list< const TVDEdge* > pointEdges;
+          edge = edge2;
+          do
+          {
+            edge = edge->next(); // Returns the CCW next edge within the cell.
+            if ( edge->is_infinite() )
+              hasInfinite = true;
+            else if ( edge->is_primary() && !inSeg.isExternal( edge ))
+              pointEdges.push_back( edge );
+          }
+          while ( edge != edge2 && !hasInfinite );
+
+          if ( hasInfinite || pointEdges.empty() )
+            break;
+          inPnt._edges.splice( inPnt._edges.end(), pointEdges );
+          inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd );
+
+          edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next();
+        }
+      } // add MA edges corresponding to concave InPoints
+
+    } // loop on inSegments to find corresponding MA edges
+
+
+    // -------------------------------------------
+    // Create Branches and BndPoints for each EDGE
+    // -------------------------------------------
+
+    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 );
+      inPoints[0]._edges.clear();
+    }
+
+    // Divide InSegment's into BndSeg's
+
+    vector< BndSeg > bndSegs;
+    bndSegs.reserve( inSegments.size() * 3 );
+
+    list< const TVDEdge* >::reverse_iterator e;
+    for ( size_t i = 0; i < inSegments.size(); ++i )
+    {
+      InSegment& inSeg = inSegments[i];
+
+      // 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
+      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:
+        vector< double > len;
+        len.push_back(0);
+        for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e )
+          len.push_back( len.back() + length( *e ));
+
+        e = inSeg._edges.rbegin();
+        for ( size_t l = 1; l < len.size(); ++e, ++l )
+        {
+          double dl = len[l] / len.back();
+          double u  = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param;
+          bndSegs.push_back( BndSeg( &inSeg, *e, u ));
+        }
+      }
+      // 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;
+    }
+
+    // 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 );
+
+
+    // Find TVDEdge's of Branches and associate them with bndSegs
+
+    vector< vector<const TVDEdge*> > branchEdges;
+    branchEdges.reserve( boundary.nbEdges() * 4 );
+
+    map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType;
+
+    int branchID = 1; // we code orientation as branchID sign
+    branchEdges.resize( branchID + 1 );
+
+    size_t i1st = 0;
+    while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
+      ++i1st;
+    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+    branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
+
+    for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
+    {
+      if ( bndSegs[i].branchID() )
+      {
+        branchID = bndSegs[i]._branchID; // with sign
+
+        if ( bndSegs[i]._branchID == -bndSegs[i-1]._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-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 ));
+      }
+      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 ( bndSegs[0]._edge )
+      {
+        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 ));
+      }
+    }
+    // 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();
+    }
+
+    // associate branchIDs and the input branch vector (arg)
+    vector< const SMESH_MAT2d::Branch* > branchByID( branchEdges.size() );
+    int nbBranches = 0;
+    for ( size_t i = 0; i < branchEdges.size(); ++i )
+    {
+      nbBranches += ( !branchEdges[i].empty() );
+    }
+    branch.resize( nbBranches );
+    for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+    {
+      if ( !branchEdges[ brID ].empty() )
+        branchByID[ brID ] = & branch[ iBr++ ];
+    }
+
+    // Fill in BndPoints of each EDGE of the boundary
+
+    size_t iSeg = 0;
+    int edgeInd = -1, dInd = 0;
+    while ( iSeg < bndSegs.size() )
+    {
+      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;
+
+      // make TVDEdge know an index of bndSegs within BndPoints
+      for ( size_t i = iSeg; i < iSegEnd; ++i )
+        if ( bndSegs[i]._edge )
+          SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge );
+
+      // parameters on EDGE
+
+      bndPoints._params.reserve( nbSegs + 1 );
+      bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param );
+
+      for ( size_t i = iSeg; i < iSegEnd; ++i )
+        bndPoints._params.push_back( bndSegs[ i ]._uLast );
+
+      // MA edges
+
+      bndPoints._maEdges.reserve( nbSegs );
+
+      for ( size_t i = iSeg; i < iSegEnd; ++i )
+      {
+        const size_t              brID = bndSegs[ i ].branchID();
+        const SMESH_MAT2d::Branch*  br = branchByID[ brID ];
+
+        if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() )
+        {
+          edgeInd += dInd;
+
+          if ( edgeInd < 0 ||
+               edgeInd >= (int) branchEdges[ brID ].size() ||
+               branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge )
+          {
+            if ( bndSegs[ i ]._branchID < 0 &&
+                 branchEdges[ brID ].back() == bndSegs[ i ]._edge )
+            {
+              edgeInd = branchEdges[ brID ].size() - 1;
+              dInd = -1;
+            }
+            else if ( bndSegs[ i ]._branchID > 0 &&
+                      branchEdges[ brID ].front() == bndSegs[ i ]._edge )
+            {
+              edgeInd = 0;
+              dInd = +1;
+            }
+            else
+            {
+              for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd )
+                if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge )
+                  break;
+              dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
+            }
+          }
+        }
+        else
+        {
+          // no MA edge, bndSeg corresponds to an end point of a branch
+          if ( bndPoints._maEdges.empty() )
+          {
+            // should not get here according to algo design
+            edgeInd = 0;
+          }
+          else
+          {
+            edgeInd = branchEdges[ brID ].size();
+            dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1;
+          }
+        }
+
+        bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd ));
+
+      } // loop on bndSegs of an EDGE
+
+      iSeg = iSegEnd;
+
+    } // loop on all bndSegs
+
+
+    // fill the branches with MA edges
+    for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID )
+      if ( !branchEdges[brID].empty() )
+      {
+        branch[ iBr ].init( branchEdges[brID], & boundary, endType );
+        iBr++;
+      }
+    // set branches to branch ends
+    for ( size_t i = 0; i < branch.size(); ++i )
+      branch[i].setBranchesToEnds( branch );
+
+    // fill branchPnt arg
+    map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end;
+    for ( size_t i = 0; i < branch.size(); ++i )
+    {
+      if ( branch[i].getEnd(0)->_branches.size() > 2 )
+        v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) ));
+      if ( branch[i].getEnd(1)->_branches.size() > 2 )
+        v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) ));
+    }
+    branchPnt.resize( v2end.size() );
+    map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin();
+    for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i )
+      branchPnt[ i ] = v2e->second;
+
+  } // makeMA()
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief MedialAxis constructor
+ *  \param [in] face - a face to create MA for
+ *  \param [in] edges - edges of the face (possibly not all) on the order they
+ *              encounter in the face boundary.
+ *  \param [in] minSegLen - minimal length of a mesh segment used to discretize
+ *              the edges. It is used to define precision of MA approximation
+ */
+//================================================================================
+
+SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face&                face,
+                                    const std::vector< TopoDS_Edge >& edges,
+                                    const double                      minSegLen,
+                                    const bool                        ignoreCorners):
+  _face( face ), _boundary( edges.size() )
+{
+  // input to construct_voronoi()
+  vector< InPoint >  inPoints;
+  vector< InSegment> inSegments;
+  if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale ))
+    return;
+
+  //inSegmentsToFile( inSegments );
+
+  // build voronoi diagram
+  construct_voronoi( inSegments.begin(), inSegments.end(), &_vd );
+
+  // make MA data
+  makeMA( _vd, inPoints, inSegments, _branch, _branchPnt, _boundary );
+}
+
+//================================================================================
+/*!
+ * \brief Return UVs of ends of MA edges of a branch
+ */
+//================================================================================
+
+void SMESH_MAT2d::MedialAxis::getPoints( const Branch&         branch,
+                                         std::vector< gp_XY >& points) const
+{
+  branch.getPoints( points, _scale );
+}
+
+//================================================================================
+/*!
+ * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE
+ *  \param [in] iGeomEdge - index of geom EDGE within a vector passed at construction
+ *  \param [in] u - parameter of the point on EDGE curve
+ *  \param [out] p - the found BranchPoint
+ *  \return bool - is OK
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge,
+                                            double            u,
+                                            BranchPoint&      p ) const
+{
+  if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
+    return false;
+
+  const BndPoints& points = _pointsPerEdge[ iEdge ];
+  const bool edgeReverse = ( points._params[0] > points._params.back() );
+
+  if ( u < ( edgeReverse ? points._params.back() : points._params[0] ))
+    u = edgeReverse ? points._params.back() : points._params[0];
+  else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) )
+    u = edgeReverse ? points._params[0] : points._params.back();
+
+  double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] );
+  int    i = int( r * double( points._maEdges.size()-1 ));
+  if ( edgeReverse )
+  {
+    while ( points._params[i  ] < u ) --i;
+    while ( points._params[i+1] > u ) ++i;
+  }
+  else
+  {
+    while ( points._params[i  ] > u ) --i;
+    while ( points._params[i+1] < u ) ++i;
+  }
+  double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] );
+
+  const std::pair< const Branch*, int >& maE = points._maEdges[ i ];
+  bool maReverse = ( maE.second < 0 );
+
+  p._branch = maE.first;
+  p._iEdge  = maE.second - 1; // countered from 1 to store sign
+  p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam;
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check if a given boundary segment is a null-length segment on a concave
+ *        boundary corner.
+ *  \param [in] iEdge - index of a geom EDGE
+ *  \param [in] iSeg - index of a boundary segment
+ *  \return bool - true if the segment is on concave corner
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const
+{
+  if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() )
+    return false;
+
+  const BndPoints& points = _pointsPerEdge[ iEdge ];
+  if ( points._params.size() >= iSeg+1 )
+    return false;
+
+  return Abs( points._params[ iEdge ] - points._params[ iEdge+1 ]) < 1e-20;
+}
+
+//================================================================================
+/*!
+ * \brief Creates a 3d curve corresponding to a Branch
+ *  \param [in] branch - the Branch
+ *  \return Adaptor3d_Curve* - the new curve the caller is to delete
+ */
+//================================================================================
+
+Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const
+{
+  Handle(Geom_Surface) surface = BRep_Tool::Surface( _face );
+  if ( surface.IsNull() )
+    return 0;
+
+  vector< gp_XY > uv;
+  branch.getPoints( uv, _scale );
+  if ( uv.size() < 2 )
+    return 0;
+
+  vector< TopoDS_Vertex > vertex( uv.size() );
+  for ( size_t i = 0; i < uv.size(); ++i )
+    vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() ));
+
+  TopoDS_Wire aWire;
+  BRep_Builder aBuilder;
+  aBuilder.MakeWire(aWire);
+  for ( size_t i = 1; i < vertex.size(); ++i )
+  {
+    TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] );
+    aBuilder.Add( aWire, edge );
+  }
+
+  // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() ))
+  //   aWire.Closed(true); // issue 0021141
+
+  return new BRepAdaptor_CompCurve( aWire );
+}
+
+//================================================================================
+/*!
+ * \brief Copy points of an EDGE
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::init( vector<const TVDEdge*>&                maEdges,
+                                const Boundary*                        boundary,
+                                map< const TVDVertex*, BranchEndType > endType )
+{
+  if ( maEdges.empty() ) return;
+
+  _boundary = boundary;
+  _maEdges.swap( maEdges );
+
+
+  _params.reserve( _maEdges.size() + 1 );
+  _params.push_back( 0. );
+  for ( size_t i = 0; i < _maEdges.size(); ++i )
+    _params.push_back( _params.back() + length( _maEdges[i] ));
+  
+  for ( size_t i = 1; i < _params.size(); ++i )
+    _params[i] /= _params.back();
+
+
+  _endPoint1._vertex = _maEdges.front()->vertex1();
+  _endPoint2._vertex = _maEdges.back ()->vertex0();
+
+  if ( endType.count( _endPoint1._vertex ))
+    _endPoint1._type = endType[ _endPoint1._vertex ];
+  if ( endType.count( _endPoint2._vertex ))
+    _endPoint2._type = endType[ _endPoint2._vertex ];
+}
+
+//================================================================================
+/*!
+ * \brief fill BranchEnd::_branches of its ends
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches )
+{
+  for ( size_t i = 0; i < branches.size(); ++i )
+  {
+    if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex ||
+         this->_endPoint1._vertex == branches[i]._endPoint2._vertex )
+      this->_endPoint1._branches.push_back( &branches[i] );
+
+    if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex ||
+         this->_endPoint2._vertex == branches[i]._endPoint2._vertex )
+      this->_endPoint2._branches.push_back( &branches[i] );
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ *  \param [in] param - [0;1] normalized param on the Branch
+ *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
+ *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
+ *  \return bool - true if the BoundaryPoint's found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(double         param,
+                                            BoundaryPoint& bp1,
+                                            BoundaryPoint& bp2 ) const
+{
+  if ( param < _params[0] || param > _params.back() )
+    return false;
+  
+  // look for an index of a MA edge by param
+  double ip = param * _params.size();
+  size_t  i = size_t( Min( int( _maEdges.size()-1), int( ip )));
+
+  while ( param < _params[i  ] ) --i;
+  while ( param > _params[i+1] ) ++i;
+
+  double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] );
+
+  return getBoundaryPoints( i, r, bp1, bp2 );
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ *  \param [in] iMAEdge - index of a MA edge within this Branch
+ *  \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge
+ *  \param [out] bp1 - BoundaryPoint on EDGE with a lower index
+ *  \param [out] bp2 - BoundaryPoint on EDGE with a higher index
+ *  \return bool - true if the BoundaryPoint's found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t    iMAEdge,
+                                            double         maEdgeParam,
+                                            BoundaryPoint& bp1,
+                                            BoundaryPoint& bp2 ) const
+{
+  if ( iMAEdge > _maEdges.size() )
+    return false;
+
+  size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] );
+  size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() );
+  size_t iSeg1  = getBndSegment( _maEdges[ iMAEdge ] );
+  size_t iSeg2  = getBndSegment( _maEdges[ iMAEdge ]->twin() );
+
+  return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) &&
+           _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 ));
+}
+
+//================================================================================
+/*!
+ * \brief Returns points on two EDGEs, equidistant from a given point of this Branch
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p,
+                                            BoundaryPoint&     bp1,
+                                            BoundaryPoint&     bp2 ) const
+{
+  return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 );
+}
+
+//================================================================================
+/*!
+ * \brief Return a parameter of a BranchPoint normalized within this Branch
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const
+{
+  if ( p._iEdge > _params.size()-1 )
+    return false;
+
+  u = ( _params[ p._iEdge   ] * ( 1 - p._edgeParam ) +
+        _params[ p._iEdge+1 ] * p._edgeParam );
+
+  return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check type of both ends
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const
+{
+  return ( _endPoint1._type == type || _endPoint2._type == type );
+}
+
+//================================================================================
+/*!
+ * \brief Returns MA points
+ *  \param [out] points - the 2d points
+ *  \param [in] scale - the scale that was used to scale the 2d space of MA
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points,
+                                     const double          scale[2]) const
+{
+  points.resize( _maEdges.size() + 1 );
+
+  points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0
+                      _maEdges[0]->vertex1()->y() / scale[1] );
+
+  for ( size_t i = 0; i < _maEdges.size(); ++i )
+    points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0],
+                          _maEdges[i]->vertex0()->y() / scale[1] );
+}
+
+//================================================================================
+/*!
+ * \brief Return indices of EDGEs equidistant from this branch
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1,
+                                        std::vector< std::size_t >& edgeIDs2 ) const
+{
+  edgeIDs1.push_back( getGeomEdge( _maEdges[0] ));
+  edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() ));
+
+  for ( size_t i = 1; i < _maEdges.size(); ++i )
+  {
+    size_t ie1 = getGeomEdge( _maEdges[i] );
+    size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
+
+    if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 );
+    if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 );
+  }
+}
+
+//================================================================================
+/*!
+ * \brief Looks for a BranchPoint position around a concave VERTEX
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >&   edgeIDs1,
+                                                   std::vector< std::size_t >&   edgeIDs2,
+                                                   std::vector< BranchPoint >&   divPoints,
+                                                   const vector<const TVDEdge*>& maEdges,
+                                                   const vector<const TVDEdge*>& maEdgesTwin,
+                                                   size_t &                    i) const
+{
+  // if there is a concave vertex between EDGEs
+  // then position of a dividing BranchPoint is undefined, it is somewhere
+  // on an arc-shaped part of the Branch around the concave vertex.
+  // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle
+  // of the arc if there is no opposite VERTEX.
+  // All null-length segments around a VERTEX belong to one of EDGEs.
+
+  BranchPoint divisionPnt;
+  divisionPnt._branch = this;
+
+  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 ] );
+  bool isConcaPrev = _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 );
+  bool isConcaNext = _boundary->IsConcaveSegment( ie1,             iSeg2 );
+  if ( !isConcaNext && !isConcaPrev )
+    return false;
+
+  bool isConcaveV = false;
+
+  int iPrev = i-1, iNext = i;
+  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() )
+    {
+      iSeg2 = getBndSegment( maEdges[ iNext ] );
+      if ( _boundary->IsConcaveSegment( ie1, iSeg2 ))
+        ++iNext;
+      else
+        break;
+    }
+    bool vertexFound = false;
+    for ( size_t iE = i+1; iE < iNext; ++iE )
+    {
+      ie2 = getGeomEdge( maEdgesTwin[iE] );
+      if ( ie2 != edgeIDs2.back() )
+      {
+        // opposite VERTEX found
+        divisionPnt._iEdge = iE;
+        divisionPnt._edgeParam = 0;
+        divPoints.push_back( divisionPnt );
+        edgeIDs1.push_back( ie1 );
+        edgeIDs2.push_back( ie2 );
+        vertexFound = true;
+        }
+    }
+    if ( vertexFound )
+    {
+      i = --iNext;
+      isConcaveV = true;
+    }
+  }
+  else if ( isConcaPrev )
+  {
+    // all null-length segments passed, find their beginning
+    while ( iPrev-1 >= 0 )
+    {
+      iSeg1 = getBndSegment( maEdges[ iPrev-1 ] );
+      if ( _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 ))
+        --iPrev;
+      else
+        break;
+    }
+  }
+
+  if ( iPrev < i-1 || iNext > i )
+  {
+    // no VERTEX on the opposite EDGE, put the Branch Point in the middle
+    double par1 = _params[ iPrev ], par2 = _params[ iNext ];
+    double midPar = 0.5 * ( par1 + par2 );
+    divisionPnt._iEdge = iPrev;
+    while ( _params[ divisionPnt._iEdge + 1 ] < midPar )
+      ++divisionPnt._iEdge;
+    divisionPnt._edgeParam =
+      ( _params[ divisionPnt._iEdge + 1 ] - midPar ) /
+      ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] );
+    divPoints.push_back( divisionPnt );
+    isConcaveV = true;
+  }
+
+  return isConcaveV;
+}
+
+//================================================================================
+/*!
+ * \brief Return indices of opposite parts of EDGEs equidistant from this branch
+ *  \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE
+ *  \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE
+ *  \param [out] divPoints - BranchPoint's located between two successive unique
+ *         pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs
+ *         of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less
+ *         than number of \a edgeIDs
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
+                                                std::vector< std::size_t >& edgeIDs2,
+                                                std::vector< BranchPoint >& divPoints) const
+{
+  edgeIDs1.clear();
+  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();
+
+  // size_t lastConcaE1 = _boundary.nbEdges();
+  // size_t lastConcaE2 = _boundary.nbEdges();
+
+  BranchPoint divisionPnt;
+  divisionPnt._branch = this;
+
+  for ( size_t i = 0; i < _maEdges.size(); ++i )
+  {
+    size_t ie1 = getGeomEdge( _maEdges[i] );
+    size_t ie2 = getGeomEdge( _maEdges[i]->twin() );
+    
+    if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
+    {
+      bool isConcaveV = false;
+      if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 )
+      {
+        isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i );
+      }
+      if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 )
+      {
+        isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i );
+      }
+
+      if ( isConcaveV )
+      {
+        ie1 = getGeomEdge( _maEdges[i] );
+        ie2 = getGeomEdge( _maEdges[i]->twin() );
+      }
+      if (( !isConcaveV ) ||
+          ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ))
+      {
+        edgeIDs1.push_back( ie1 );
+        edgeIDs2.push_back( ie2 );
+      }
+      if ( divPoints.size() < edgeIDs1.size() - 1 )
+      {
+        divisionPnt._iEdge = i;
+        divisionPnt._edgeParam = 0;
+        divPoints.push_back( divisionPnt );
+      }
+
+    } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )
+  } // loop on _maEdges
+}
+
+//================================================================================
+/*!
+ * \brief Store data of boundary segments in TVDEdge
+ */
+//================================================================================
+
+void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge )
+{
+  if ( maEdge ) maEdge->cell()->color( geomIndex );
+}
+std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge )
+{
+  return maEdge ? maEdge->cell()->color() : std::string::npos;
+}
+void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge )
+{
+  if ( maEdge ) maEdge->color( segIndex );
+}
+std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge )
+{
+  return maEdge ? maEdge->color() : std::string::npos;
+}
+
+//================================================================================
+/*!
+ * \brief Returns a boundary point on a given EDGE
+ *  \param [in] iEdge - index of the EDGE within MedialAxis
+ *  \param [in] iSeg - index of a boundary segment within this Branch
+ *  \param [in] u - [0;1] normalized param within \a iSeg-th segment
+ *  \param [out] bp - the found BoundaryPoint
+ *  \return bool - true if the BoundaryPoint is found
+ */
+//================================================================================
+
+bool SMESH_MAT2d::Boundary::getPoint( std::size_t    iEdge,
+                                      std::size_t    iSeg,
+                                      double         u,
+                                      BoundaryPoint& bp ) const
+{
+  if ( iEdge >= _pointsPerEdge.size() )
+    return false;
+  if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() )
+    return false;
+
+  // This method is called by Branch that can have an opposite orientation,
+  // hence u is inverted depending on orientation coded as a sign of _maEdge index
+  bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 );
+  if ( isReverse )
+    u = 1. - u;
+
+  double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ];
+  double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ];
+
+  bp._param = p0 * ( 1. - u ) + p1 * u;
+  bp._edgeIndex = iEdge;
+
+  return true;
+}
+
diff --git a/src/SMESHUtils/SMESH_MAT2d.hxx b/src/SMESHUtils/SMESH_MAT2d.hxx
new file mode 100644
index 000000000..7e1061d05
--- /dev/null
+++ b/src/SMESHUtils/SMESH_MAT2d.hxx
@@ -0,0 +1,224 @@
+// Copyright (C) 2007-2015  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
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : SMESH_MAT2d.hxx
+// Created   : Thu May 28 17:49:53 2015
+// Author    : Edward AGAPOV (eap)
+
+#ifndef __SMESH_MAT2d_HXX__
+#define __SMESH_MAT2d_HXX__
+
+#include "SMESH_Utils.hxx"
+
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Edge.hxx>
+
+#include <vector>
+#include <map>
+
+#include <boost/polygon/polygon.hpp>
+#include <boost/polygon/voronoi.hpp>
+
+class Adaptor3d_Curve;
+
+// Medial Axis Transform 2D
+namespace SMESH_MAT2d
+{
+  class MedialAxis; // MedialAxis is the entry point
+  class Branch;
+  class BranchEnd;
+  class Boundary;
+  struct BoundaryPoint;
+
+  typedef boost::polygon::voronoi_diagram<double> TVD;
+  typedef TVD::cell_type                          TVDCell;
+  typedef TVD::edge_type                          TVDEdge;
+  typedef TVD::vertex_type                        TVDVertex;
+
+  //-------------------------------------------------------------------------------------
+  // type of Branch end point
+  enum BranchEndType { BE_UNDEF,
+                       BE_ON_VERTEX, // branch ends at a convex VRTEX
+                       BE_BRANCH_POINT, // branch meats 2 or more other branches
+                       BE_END // branch end equidistant from several adjacent segments
+  };
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief End point of MA Branch
+   */
+  struct SMESHUtils_EXPORT BranchEnd
+  {
+    const TVDVertex*             _vertex;
+    BranchEndType                _type;
+    std::vector< const Branch* > _branches;
+
+    BranchEnd(): _vertex(0), _type( BE_UNDEF ) {}
+  };
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Point on MA Branch
+   */
+  struct SMESHUtils_EXPORT BranchPoint
+  {
+    const Branch* _branch;
+    std::size_t   _iEdge; // MA edge index within the branch
+    double        _edgeParam; // normalized param within the MA edge
+
+    BranchPoint(): _branch(0), _iEdge(0), _edgeParam(-1) {}
+  };
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Branch is a set of MA edges enclosed between branch points and/or MA ends.
+   *        It's main feature is to return two BoundaryPoint's per a point on it.
+   */
+  class SMESHUtils_EXPORT Branch
+  {
+  public:
+    bool getBoundaryPoints(double param, BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+    bool getBoundaryPoints(std::size_t iMAEdge, double maEdgeParam,
+                           BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+    bool getBoundaryPoints(const BranchPoint& p,
+                           BoundaryPoint& bp1, BoundaryPoint& bp2 ) const;
+    bool getParameter(const BranchPoint& p, double & u ) const;
+
+    std::size_t      nbEdges() const { return _maEdges.size(); }
+
+    const BranchEnd* getEnd(bool the2nd) const { return & ( the2nd ? _endPoint2 : _endPoint1 ); }
+
+    bool hasEndOfType(BranchEndType type) const;
+
+    void getPoints( std::vector< gp_XY >& points, const double scale[2]) const;
+
+    void getGeomEdges( std::vector< std::size_t >& edgeIDs1,
+                       std::vector< std::size_t >& edgeIDs2 ) const;
+
+    void getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1,
+                               std::vector< std::size_t >& edgeIDs2,
+                               std::vector< BranchPoint >& divPoints) const;
+
+    // construction
+    void init( std::vector<const TVDEdge*>&                maEdges,
+               const Boundary*                             boundary,
+               std::map< const TVDVertex*, BranchEndType > endType);
+    void setBranchesToEnds( const std::vector< Branch >&   branches);
+
+    static void setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge );
+    static std::size_t getGeomEdge( const TVDEdge* maEdge );
+    static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge );
+    static std::size_t getBndSegment( const TVDEdge* maEdge );
+
+  private:
+
+    bool addDivPntForConcaVertex( std::vector< std::size_t >&        edgeIDs1,
+                                  std::vector< std::size_t >&        edgeIDs2,
+                                  std::vector< BranchPoint >&        divPoints,
+                                  const std::vector<const TVDEdge*>& maEdges,
+                                  const std::vector<const TVDEdge*>& maEdgesTwin,
+                                  size_t &                           i) const;
+
+    // association of _maEdges with boundary segments is stored in this way:
+    // index of an EDGE:           TVDEdge->cell()->color()
+    // index of a segment on EDGE: TVDEdge->color()
+    std::vector<const TVDEdge*> _maEdges; // MA edges ending at points located at _params
+    std::vector<double>  _params; // params of points on MA, normalized [0;1] within this branch
+    const Boundary*      _boundary; // face boundary
+    BranchEnd            _endPoint1;
+    BranchEnd            _endPoint2;
+  };
+
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Data of a discretized EDGE allowing to get a point on MA by a parameter on EDGE
+   */
+  struct BndPoints
+  {
+    std::vector< double > _params; // params of discretization points on an EDGE
+    std::vector< std::pair< const Branch*, int > > _maEdges; /* index of TVDEdge in branch;
+                                                                index sign means orientation;
+                                                                index == Branch->nbEdges() means
+                                                                end point of a Branch */
+  };
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Face boundary is discretized so that each its segment to correspond to
+   *        an edge of MA
+   */
+  class Boundary
+  {
+  public:
+
+    Boundary( std::size_t nbEdges ): _pointsPerEdge( nbEdges ) {}
+    BndPoints&  getPoints( std::size_t iEdge ) { return _pointsPerEdge[ iEdge ]; }
+    std::size_t nbEdges() const { return _pointsPerEdge.size(); }
+
+    bool getPoint( std::size_t iEdge, std::size_t iSeg, double u, BoundaryPoint& bp ) const;
+
+    bool getBranchPoint( const std::size_t iEdge, double u, BranchPoint& p ) const;
+
+    bool IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const;
+
+  private:
+    std::vector< BndPoints > _pointsPerEdge;
+  };
+
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Point on FACE boundary
+   */
+  struct SMESHUtils_EXPORT BoundaryPoint
+  {
+    std::size_t _edgeIndex; // index of an EDGE in a sequence passed to MedialAxis()
+    double      _param;     // parameter of this EDGE
+  };
+  //-------------------------------------------------------------------------------------
+  /*!
+   * \brief Medial axis (MA) is defined as the loci of centres of locally
+   *        maximal balls inside 2D representation of a face. This class
+   *        implements a piecewise approximation of MA.
+   */
+  class SMESHUtils_EXPORT MedialAxis
+  {
+  public:
+    MedialAxis(const TopoDS_Face&                face,
+               const std::vector< TopoDS_Edge >& edges,
+               const double                      minSegLen,
+               const bool                        ignoreCorners = false );
+    const Boundary&                        getBoundary() const { return _boundary; }
+    const std::vector< Branch >&           getBranches() const { return _branch; }
+    const std::vector< const BranchEnd* >& getBranchPoints() const { return _branchPnt; }
+
+    void getPoints( const Branch& branch, std::vector< gp_XY >& points) const;
+    Adaptor3d_Curve* make3DCurve(const Branch& branch) const;
+
+  private:
+
+  private:
+    TopoDS_Face                     _face;
+    TVD                             _vd;
+    std::vector< Branch >           _branch;
+    std::vector< const BranchEnd* > _branchPnt;
+    Boundary                        _boundary;
+    double                          _scale[2];
+  };
+
+}
+
+#endif
diff --git a/src/SMESH_SWIG/StdMeshersBuilder.py b/src/SMESH_SWIG/StdMeshersBuilder.py
index 41e0c2269..3a07810ad 100644
--- a/src/SMESH_SWIG/StdMeshersBuilder.py
+++ b/src/SMESH_SWIG/StdMeshersBuilder.py
@@ -42,6 +42,8 @@ Hexa        = "Hexa_3D"
 QUADRANGLE  = "Quadrangle_2D"
 ## Algorithm type: Radial Quadrangle 1D-2D algorithm, see StdMeshersBuilder_RadialQuadrangle1D2D
 RADIAL_QUAD = "RadialQuadrangle_1D2D"
+## Algorithm type: Quadrangle (Medial Axis Projection) 1D-2D algorithm, see StdMeshersBuilder_QuadMA_1D2D
+QUAD_MA_PROJ = "QuadFromMedialAxis_1D2D"
 
 # import items of enums
 for e in StdMeshers.QuadType._items: exec('%s = StdMeshers.%s'%(e,e))
@@ -455,9 +457,6 @@ class StdMeshersBuilder_Segment_Python(Mesh_Algorithm):
     algoType   = PYTHON
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates tetrahedron 3D algorithm for solids"
-    ## doc string of the method
-    #  @internal
     docHelper  = "Creates segment 1D algorithm for edges"
 
     ## Private constructor.
@@ -856,7 +855,7 @@ class StdMeshersBuilder_Projection1D2D(StdMeshersBuilder_Projection2D):
     algoType   = "Projection_1D2D"
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates projection 1D-2D algorithm for edges and faces"
+    docHelper  = "Creates projection 1D-2D algorithm for faces"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1115,7 +1114,7 @@ class StdMeshersBuilder_RadialPrism3D(StdMeshersBuilder_Prism3D):
     algoType   = "RadialPrism_3D"
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates prism 3D algorithm for volumes"
+    docHelper  = "Creates Raial Prism 3D algorithm for volumes"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1147,7 +1146,7 @@ class StdMeshersBuilder_RadialQuadrangle1D2D(Mesh_Algorithm):
     algoType   = RADIAL_QUAD
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates quadrangle 1D-2D algorithm for triangular faces"
+    docHelper  = "Creates quadrangle 1D-2D algorithm for faces having a shape of disk or a disk segment"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1258,6 +1257,35 @@ class StdMeshersBuilder_RadialQuadrangle1D2D(Mesh_Algorithm):
 
     pass # end of StdMeshersBuilder_RadialQuadrangle1D2D class
 
+## Defines a Quadrangle (Medial Axis Projection) 1D-2D algorithm
+# 
+#  It is created by calling smeshBuilder.Mesh.Quadrangle(smeshBuilder.QUAD_MA_PROJ,geom=0)
+#
+#  @ingroup l2_algos_quad_ma
+class StdMeshersBuilder_QuadMA_1D2D(Mesh_Algorithm):
+
+    ## name of the dynamic method in smeshBuilder.Mesh class
+    #  @internal
+    meshMethod = "Quadrangle"
+    ## type of algorithm used with helper function in smeshBuilder.Mesh class
+    #  @internal
+    algoType   = QUAD_MA_PROJ
+    ## doc string of the method
+    #  @internal
+    docHelper  = "Creates quadrangle 1D-2D algorithm for faces"
+
+    ## Private constructor.
+    #  @param mesh parent mesh object algorithm is assigned to
+    #  @param geom geometry (shape/sub-shape) algorithm is assigned to;
+    #              if it is @c 0 (default), the algorithm is assigned to the main shape
+    def __init__(self, mesh, geom=0):
+        Mesh_Algorithm.__init__(self)
+        self.Create(mesh, geom, self.algoType)
+        pass
+
+    pass
+
+
 ## Defines a Use Existing Elements 1D algorithm
 #
 #  It is created by calling smeshBuilder.Mesh.UseExisting1DElements(geom=0)
@@ -1327,7 +1355,7 @@ class StdMeshersBuilder_UseExistingElements_1D2D(Mesh_Algorithm):
     isDefault  = True
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates 1D-2D algorithm for edges/faces with reusing of existing mesh elements"
+    docHelper  = "Creates 1D-2D algorithm for faces with reusing of existing mesh elements"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1375,7 +1403,7 @@ class StdMeshersBuilder_Cartesian_3D(Mesh_Algorithm):
     isDefault  = True
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates body fitting 3D algorithm for volumes"
+    docHelper  = "Creates Body Fitting 3D algorithm for volumes"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1505,7 +1533,7 @@ class StdMeshersBuilder_UseExisting_1D(Mesh_Algorithm):
     algoType   = "UseExisting_1D"
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates 1D algorithm for edges with reusing of existing mesh elements"
+    docHelper  = "Creates 1D algorithm allowing batch meshing of edges"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
@@ -1533,7 +1561,7 @@ class StdMeshersBuilder_UseExisting_2D(Mesh_Algorithm):
     algoType   = "UseExisting_2D"
     ## doc string of the method
     #  @internal
-    docHelper  = "Creates 2D algorithm for faces with reusing of existing mesh elements"
+    docHelper  = "Creates 2D algorithm allowing batch meshing of faces"
 
     ## Private constructor.
     #  @param mesh parent mesh object algorithm is assigned to
diff --git a/src/StdMeshers/CMakeLists.txt b/src/StdMeshers/CMakeLists.txt
index 46b756b75..8714c1b4f 100644
--- a/src/StdMeshers/CMakeLists.txt
+++ b/src/StdMeshers/CMakeLists.txt
@@ -130,6 +130,7 @@ SET(StdMeshers_HEADERS
   StdMeshers_Projection_1D2D.hxx
   StdMeshers_CartesianParameters3D.hxx
   StdMeshers_Cartesian_3D.hxx
+  StdMeshers_QuadFromMedialAxis_1D2D.hxx
 )
 
 IF(SALOME_SMESH_ENABLE_MEFISTO)
@@ -193,6 +194,7 @@ SET(StdMeshers_SOURCES
   StdMeshers_CartesianParameters3D.cxx
   StdMeshers_Cartesian_3D.cxx
   StdMeshers_Adaptive1D.cxx
+  StdMeshers_QuadFromMedialAxis_1D2D.cxx
 )
 
 IF(SALOME_SMESH_ENABLE_MEFISTO)
diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx
index 0616b630d..264189170 100644
--- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx
+++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx
@@ -696,7 +696,7 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh &        aMesh,
   double xmax = -1.e300;
   double ymin = 1.e300;
   double ymax = -1.e300;
-  int nbp = 23;
+  const int nbp = 23;
   scalex = 1;
   scaley = 1;
 
@@ -720,13 +720,8 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh &        aMesh,
         ymin = p.Y();
       if (p.Y() > ymax)
         ymax = p.Y();
-      //    MESSAGE(" "<< f<<" "<<l<<" "<<param<<" "<<xmin<<" "<<xmax<<" "<<ymin<<" "<<ymax);
     }
   }
-  //   SCRUTE(xmin);
-  //   SCRUTE(xmax);
-  //   SCRUTE(ymin);
-  //   SCRUTE(ymax);
   double xmoy = (xmax + xmin) / 2.;
   double ymoy = (ymax + ymin) / 2.;
   double xsize = xmax - xmin;
@@ -754,23 +749,14 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh &        aMesh,
   }
   scalex = length_x / xsize;
   scaley = length_y / ysize;
-//   SCRUTE(xsize);
-//   SCRUTE(ysize);
   double xyratio = xsize*scalex/(ysize*scaley);
   const double maxratio = 1.e2;
-  //SCRUTE(xyratio);
   if (xyratio > maxratio) {
-    SCRUTE( scaley );
     scaley *= xyratio / maxratio;
-    SCRUTE( scaley );
   }
   else if (xyratio < 1./maxratio) {
-    SCRUTE( scalex );
     scalex *= 1 / xyratio / maxratio;
-    SCRUTE( scalex );
   }
-  ASSERT(scalex);
-  ASSERT(scaley);
 }
 
 // namespace
diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx
new file mode 100644
index 000000000..a6d372d5f
--- /dev/null
+++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx
@@ -0,0 +1,1251 @@
+// Copyright (C) 2007-2015  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
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : StdMeshers_QuadFromMedialAxis_1D2D.cxx
+// Created   : Wed Jun  3 17:33:45 2015
+// Author    : Edward AGAPOV (eap)
+
+#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx"
+
+#include "SMESH_Block.hxx"
+#include "SMESH_Gen.hxx"
+#include "SMESH_MAT2d.hxx"
+#include "SMESH_Mesh.hxx"
+#include "SMESH_MesherHelper.hxx"
+#include "SMESH_ProxyMesh.hxx"
+#include "SMESH_subMesh.hxx"
+#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_Regular_1D.hxx"
+#include "StdMeshers_ViscousLayers2D.hxx"
+
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepTools.hxx>
+#include <BRep_Tool.hxx>
+#include <GeomAPI_Interpolate.hxx>
+#include <Geom_Surface.hxx>
+#include <Precision.hxx>
+#include <TColgp_HArray1OfPnt.hxx>
+#include <TopExp.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <gp_Pnt.hxx>
+
+#include <list>
+#include <vector>
+
+//================================================================================
+/*!
+ * \brief 1D algo
+ */
+class StdMeshers_QuadFromMedialAxis_1D2D::Algo1D : public StdMeshers_Regular_1D
+{
+public:
+  Algo1D(int studyId, SMESH_Gen* gen):
+    StdMeshers_Regular_1D( gen->GetANewId(), studyId, gen )
+  {
+  }
+  void SetSegmentLength( double len )
+  {
+    _value[ BEG_LENGTH_IND ] = len;
+    _value[ PRECISION_IND ] = 1e-7;
+    _hypType = LOCAL_LENGTH;
+  }
+};
+ 
+//================================================================================
+/*!
+ * \brief Constructor sets algo features
+ */
+//================================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int        hypId,
+                                                                       int        studyId,
+                                                                       SMESH_Gen* gen)
+  : StdMeshers_Quadrangle_2D(hypId, studyId, gen),
+    _regular1D( 0 )
+{
+  _name = "QuadFromMedialAxis_1D2D";
+  _shapeType = (1 << TopAbs_FACE);
+  _onlyUnaryInput          = true;  // FACE by FACE so far
+  _requireDiscreteBoundary = false; // make 1D by myself
+  _supportSubmeshes        = true; // make 1D by myself
+  _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
+  _neededLowerHyps[ 2 ]    = true;  // suppress warning on hiding a global 2D algo
+  _compatibleHypothesis.clear();
+  //_compatibleHypothesis.push_back("ViscousLayers2D");
+}
+
+//================================================================================
+/*!
+ * \brief Destructor
+ */
+//================================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D::~StdMeshers_QuadFromMedialAxis_1D2D()
+{
+  delete _regular1D;
+  _regular1D = 0;
+}
+
+//================================================================================
+/*!
+ * \brief Check if needed hypotheses are present
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh&         aMesh,
+                                                         const TopoDS_Shape& aShape,
+                                                         Hypothesis_Status&  aStatus)
+{
+  return true; // does not require hypothesis
+}
+
+namespace
+{
+  //================================================================================
+  /*!
+   * \brief Temporary mesh
+   */
+  struct TmpMesh : public SMESH_Mesh
+  {
+    TmpMesh()
+    {
+      _myMeshDS = new SMESHDS_Mesh(/*id=*/0, /*isEmbeddedMode=*/true);
+    }
+  };
+
+  //================================================================================
+  /*!
+   * \brief Select two EDGEs from a map, either mapped to least values or to max values
+   */
+  //================================================================================
+
+  // template< class TVal2EdgesMap >
+  // void getTwo( bool                 least,
+  //              TVal2EdgesMap&       map,
+  //              vector<TopoDS_Edge>& twoEdges,
+  //              vector<TopoDS_Edge>& otherEdges)
+  // {
+  //   twoEdges.clear();
+  //   otherEdges.clear();
+  //   if ( least )
+  //   {
+  //     TVal2EdgesMap::iterator i = map.begin();
+  //     twoEdges.push_back( i->second );
+  //     twoEdges.push_back( ++i->second );
+  //     for ( ; i != map.end(); ++i )
+  //       otherEdges.push_back( i->second );
+  //   }
+  //   else
+  //   {
+  //     TVal2EdgesMap::reverse_iterator i = map.rbegin();
+  //     twoEdges.push_back( i->second );
+  //     twoEdges.push_back( ++i->second );
+  //     for ( ; i != map.rend(); ++i )
+  //       otherEdges.push_back( i->second );
+  //   }
+  //   TopoDS_Vertex v;
+  //   if ( TopExp::CommonVertex( twoEdges[0], twoEdges[1], v ))
+  //   {
+  //     twoEdges.clear(); // two EDGEs must not be connected
+  //     otherEdges.clear();
+  //   }
+  // }
+
+  //================================================================================
+  /*!
+   * \brief Finds out a minimal segment length given EDGEs will be divided into.
+   *        This length is further used to discretize the Medial Axis
+   */
+  //================================================================================
+
+  double getMinSegLen(SMESH_MesherHelper&        theHelper,
+                      const vector<TopoDS_Edge>& theEdges)
+  {
+    TmpMesh tmpMesh;
+    SMESH_Mesh* mesh = theHelper.GetMesh();
+
+    vector< SMESH_Algo* > algos( theEdges.size() );
+    for ( size_t i = 0; i < theEdges.size(); ++i )
+    {
+      SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
+      algos[i] = sm->GetAlgo();
+    }
+
+    const int nbSegDflt = mesh->GetGen()->GetDefaultNbSegments();
+    double minSegLen    = Precision::Infinite();
+
+    for ( size_t i = 0; i < theEdges.size(); ++i )
+    {
+      SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] );
+      if ( SMESH_Algo::IsStraight( theEdges[i], /*degenResult=*/true ))
+        continue;
+      // get algo
+      size_t iOpp = ( theEdges.size() == 4 ? (i+2)%4 : i );
+      SMESH_Algo*  algo = sm->GetAlgo();
+      if ( !algo ) algo = algos[ iOpp ];
+      // get hypo
+      SMESH_Hypothesis::Hypothesis_Status status = SMESH_Hypothesis::HYP_MISSING;
+      if ( algo )
+      {
+        if ( !algo->CheckHypothesis( *mesh, theEdges[i], status ))
+          algo->CheckHypothesis( *mesh, theEdges[iOpp], status );
+      }
+      // compute
+      if ( status != SMESH_Hypothesis::HYP_OK )
+      {
+        minSegLen = Min( minSegLen, SMESH_Algo::EdgeLength( theEdges[i] ) / nbSegDflt );
+      }
+      else
+      {
+        tmpMesh.Clear();
+        tmpMesh.ShapeToMesh( TopoDS_Shape());
+        tmpMesh.ShapeToMesh( theEdges[i] );
+        try {
+          mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes
+          if ( !algo->Compute( tmpMesh, theEdges[i] ))
+            continue;
+        }
+        catch (...) {
+          continue;
+        }
+        SMDS_EdgeIteratorPtr segIt = tmpMesh.GetMeshDS()->edgesIterator();
+        while ( segIt->more() )
+        {
+          const SMDS_MeshElement* seg = segIt->next();
+          double len = SMESH_TNodeXYZ( seg->GetNode(0) ).Distance( seg->GetNode(1) );
+          minSegLen = Min( minSegLen, len );
+        }
+      }
+    }
+    if ( Precision::IsInfinite( minSegLen ))
+      minSegLen = mesh->GetShapeDiagonalSize() / nbSegDflt;
+
+    return minSegLen;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Returns EDGEs located between two VERTEXes at which given MA branches end
+   *  \param [in] br1 - one MA branch
+   *  \param [in] br2 - one more MA branch
+   *  \param [in] allEdges - all EDGEs of a FACE
+   *  \param [out] shortEdges - the found EDGEs
+   *  \return bool - is OK or not
+   */
+  //================================================================================
+
+  bool getConnectedEdges( const SMESH_MAT2d::Branch* br1,
+                          const SMESH_MAT2d::Branch* br2,
+                          const vector<TopoDS_Edge>& allEdges,
+                          vector<TopoDS_Edge>&       shortEdges)
+  {
+    vector< size_t > edgeIDs[4];
+    br1->getGeomEdges( edgeIDs[0], edgeIDs[1] );
+    br2->getGeomEdges( edgeIDs[2], edgeIDs[3] );
+
+    // EDGEs returned by a Branch form a connected chain with a VERTEX where
+    // the Branch ends at the chain middle. One of end EDGEs of the chain is common
+    // with either end EDGE of the chain of the other Branch, or the chains are connected
+    // at a common VERTEX;
+
+    // Get indices of end EDGEs of the branches
+    bool vAtStart1 = ( br1->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
+    bool vAtStart2 = ( br2->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX );
+    size_t iEnd[4] = {
+      vAtStart1 ? edgeIDs[0].back() : edgeIDs[0][0],
+      vAtStart1 ? edgeIDs[1].back() : edgeIDs[1][0],
+      vAtStart2 ? edgeIDs[2].back() : edgeIDs[2][0],
+      vAtStart2 ? edgeIDs[3].back() : edgeIDs[3][0]
+    };
+
+    set< size_t > connectedIDs;
+    TopoDS_Vertex vCommon;
+    // look for the same EDGEs
+    for ( int i = 0; i < 2; ++i )
+      for ( int j = 2; j < 4; ++j )
+        if ( iEnd[i] == iEnd[j] )
+        {
+          connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
+          connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
+          i = j = 4;
+        }
+    if ( connectedIDs.empty() )
+      // look for connected EDGEs
+      for ( int i = 0; i < 2; ++i )
+        for ( int j = 2; j < 4; ++j )
+          if ( TopExp::CommonVertex( allEdges[ iEnd[i]], allEdges[ iEnd[j]], vCommon ))
+          {
+            connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() );
+            connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() );
+            i = j = 4;
+          }
+    if ( connectedIDs.empty() ||                     // nothing
+         allEdges.size() - connectedIDs.size() < 2 ) // too many
+      return false;
+
+    // set shortEdges in the order as in allEdges
+    if ( connectedIDs.count( 0 ) &&
+         connectedIDs.count( allEdges.size()-1 ))
+    {
+      size_t iE = allEdges.size()-1;
+      while ( connectedIDs.count( iE-1 ))
+        --iE;
+      for ( size_t i = 0; i < connectedIDs.size(); ++i )
+      {
+        shortEdges.push_back( allEdges[ iE ]);
+        iE = ( iE + 1 ) % allEdges.size();
+      }
+    }
+    else
+    {
+      set< size_t >::iterator i = connectedIDs.begin();
+      for ( ; i != connectedIDs.end(); ++i )
+        shortEdges.push_back( allEdges[ *i ]);
+    }
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Find EDGEs to discretize using projection from MA
+   *  \param [in] theFace - the FACE to be meshed
+   *  \param [in] theWire - ordered EDGEs of the FACE
+   *  \param [out] theSinuEdges - the EDGEs to discretize using projection from MA
+   *  \param [out] theShortEdges - other EDGEs
+   *  \return bool - OK or not
+   *
+   * Is separate all EDGEs into four sides of a quadrangle connected in the order:
+   * theSinuEdges[0], theShortEdges[0], theSinuEdges[1], theShortEdges[1]
+   */
+  //================================================================================
+
+  bool getSinuousEdges( SMESH_MesherHelper& theHelper,
+                        const TopoDS_Face&  theFace,
+                        list<TopoDS_Edge>&  theWire,
+                        vector<TopoDS_Edge> theSinuEdges[2],
+                        vector<TopoDS_Edge> theShortEdges[2])
+  {
+    theSinuEdges[0].clear();
+    theSinuEdges[1].clear();
+    theShortEdges[0].clear();
+    theShortEdges[1].clear();
+
+    vector<TopoDS_Edge> allEdges( theWire.begin(), theWire.end() );
+    const size_t nbEdges = allEdges.size();
+    if ( nbEdges < 4 )
+      return false;
+
+    // create MedialAxis to find short edges by analyzing MA branches
+    double minSegLen = getMinSegLen( theHelper, allEdges );
+    SMESH_MAT2d::MedialAxis ma( theFace, allEdges, minSegLen );
+
+    // in an initial request case, theFace represents a part of a river with almost parallel banks
+    // so there should be two branch points
+    using SMESH_MAT2d::BranchEnd;
+    using SMESH_MAT2d::Branch;
+    const vector< const BranchEnd* >& braPoints = ma.getBranchPoints();
+    if ( braPoints.size() < 2 )
+      return false;
+    TopTools_MapOfShape shortMap;
+    size_t nbBranchPoints = 0;
+    for ( size_t i = 0; i < braPoints.size(); ++i )
+    {
+      vector< const Branch* > vertBranches; // branches with an end on VERTEX
+      for ( size_t ib = 0; ib < braPoints[i]->_branches.size(); ++ib )
+      {
+        const Branch* branch = braPoints[i]->_branches[ ib ];
+        if ( branch->hasEndOfType( SMESH_MAT2d::BE_ON_VERTEX ))
+          vertBranches.push_back( branch );
+      }
+      if ( vertBranches.size() != 2 || braPoints[i]->_branches.size() != 3)
+        continue;
+
+      // get common EDGEs of two branches
+      if ( !getConnectedEdges( vertBranches[0], vertBranches[1],
+                               allEdges, theShortEdges[ nbBranchPoints > 0 ] ))
+        return false;
+
+      for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS )
+        shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]);
+
+      ++nbBranchPoints;
+    }
+
+    if ( nbBranchPoints != 2 )
+      return false;
+
+    // add to theSinuEdges all edges that are not theShortEdges
+    vector< vector<TopoDS_Edge> > sinuEdges(1);
+    TopoDS_Vertex vCommon;
+    for ( size_t i = 0; i < allEdges.size(); ++i )
+    {
+      if ( !shortMap.Contains( allEdges[i] ))
+      {
+        if ( !sinuEdges.back().empty() )
+          if ( !TopExp::CommonVertex( sinuEdges.back().back(), allEdges[ i ], vCommon ))
+            sinuEdges.resize( sinuEdges.size() + 1 );
+
+        sinuEdges.back().push_back( allEdges[i] );
+      }
+    }
+    if ( sinuEdges.size() == 3 )
+    {
+      if ( !TopExp::CommonVertex( sinuEdges.back().back(), sinuEdges[0][0], vCommon ))
+        return false;
+      vector<TopoDS_Edge>& last = sinuEdges.back();
+      last.insert( last.end(), sinuEdges[0].begin(), sinuEdges[0].end() );
+      sinuEdges[0].swap( last );
+      sinuEdges.resize( 2 );
+    }
+    if ( sinuEdges.size() != 2 )
+      return false;
+
+    theSinuEdges[0].swap( sinuEdges[0] );
+    theSinuEdges[1].swap( sinuEdges[1] );
+
+    if ( !TopExp::CommonVertex( theSinuEdges[0].back(), theShortEdges[0][0], vCommon ) ||
+         !vCommon.IsSame( theHelper.IthVertex( 1, theSinuEdges[0].back() )))
+      theShortEdges[0].swap( theShortEdges[1] );
+
+    return ( theShortEdges[0].size() > 0 && theShortEdges[1].size() > 0 &&
+             theSinuEdges [0].size() > 0 && theSinuEdges [1].size() > 0 );
+
+    // the sinuous EDGEs can be composite and C0 continuous,
+    // therefor we use a complex criterion to find TWO short non-sinuous EDGEs
+    // and the rest EDGEs will be treated as sinuous.
+    // A short edge should have the following features:
+  // a) straight
+  // b) short
+  // c) with convex corners at ends
+  // d) far from the other short EDGE
+
+  // vector< double > isStraightEdge( nbEdges, 0 ); // criterion value
+
+  // // a0) evaluate continuity
+  // const double contiWgt = 0.5; // weight of continuity in the criterion
+  // multimap< int, TopoDS_Edge > continuity;
+  // for ( size_t i = 0; i < nbEdges; ++I )
+    // {
+    //   BRepAdaptor_Curve curve( allEdges[i] );
+    //   GeomAbs_Shape C = GeomAbs_CN;
+    //   try:
+    //     C = curve.Continuity(); // C0, G1, C1, G2, C2, C3, CN
+    //   catch ( Standard_Failure ) {}
+    //   continuity.insert( make_pair( C, allEdges[i] ));
+    //   isStraight[i] += double( C ) / double( CN ) * contiWgt;
+    // }
+
+    // // try to choose by continuity
+    // int mostStraight = (int) continuity.rbegin()->first;
+    // int lessStraight = (int) continuity.begin()->first;
+    // if ( mostStraight != lessStraight )
+    // {
+    //   int nbStraight = continuity.count( mostStraight );
+    //   if ( nbStraight == 2 )
+    //   {
+    //     getTwo( /*least=*/false, continuity, theShortEdges, theSinuEdges );
+    //   }
+    //   else if ( nbStraight == 3 && nbEdges == 4 )
+    //   {
+    //     theSinuEdges.push_back( continuity.begin()->second );
+    //     vector<TopoDS_Edge>::iterator it =
+    //       std::find( allEdges.begin(), allEdges.end(), theSinuEdges[0] );
+    //     int i = std::distance( allEdges.begin(), it );
+    //     theSinuEdges .push_back( allEdges[( i+2 )%4 ]);
+    //     theShortEdges.push_back( allEdges[( i+1 )%4 ]);
+    //     theShortEdges.push_back( allEdges[( i+3 )%4 ]);
+    //   }
+    //   if ( theShortEdges.size() == 2 )
+    //     return true;
+    // }
+
+    // // a) curvature; evaluate aspect ratio
+    // {
+    //   const double curvWgt = 0.5;
+    //   for ( size_t i = 0; i < nbEdges; ++I )
+    //   {
+    //     BRepAdaptor_Curve curve( allEdges[i] );
+    //     double curvature = 1;
+    //     if ( !curve.IsClosed() )
+    //     {
+    //       const double f = curve.FirstParameter(), l = curve.LastParameter();
+    //       gp_Pnt pf = curve.Value( f ), pl = curve.Value( l );
+    //       gp_Lin line( pf, pl.XYZ() - pf.XYZ() );
+    //       double distMax = 0;
+    //       for ( double u = f; u < l; u += (l-f)/30. )
+    //         distMax = Max( distMax, line.SquareDistance( curve.Value( u )));
+    //       curvature = Sqrt( distMax ) / ( pf.Distance( pl ));
+    //     }
+    //     isStraight[i] += curvWgt / (              curvature + 1e-20 );
+    //   }
+    // }
+    // // b) length
+    // {
+    //   const double lenWgt = 0.5;
+    //   for ( size_t i = 0; i < nbEdges; ++I )
+    //   {
+    //     double length = SMESH_Algo::Length( allEdges[i] );
+    //     if ( length > 0 )
+    //       isStraight[i] += lenWgt / length;
+    //   }
+    // }
+    // // c) with convex corners at ends
+    // {
+    //   const double cornerWgt = 0.25;
+    //   for ( size_t i = 0; i < nbEdges; ++I )
+    //   {
+    //     double convex = 0;
+    //     int iPrev = SMESH_MesherHelper::WrapIndex( int(i)-1, nbEdges );
+    //     int iNext = SMESH_MesherHelper::WrapIndex( int(i)+1, nbEdges );
+    //     TopoDS_Vertex v = helper.IthVertex( 0, allEdges[i] );
+    //     double angle = SMESH_MesherHelper::GetAngle( allEdges[iPrev], allEdges[i], theFace, v );
+    //     if ( angle < M_PI ) // [-PI; PI]
+    //       convex += ( angle + M_PI ) / M_PI / M_PI;
+    //     v = helper.IthVertex( 1, allEdges[i] );
+    //     angle = SMESH_MesherHelper::GetAngle( allEdges[iNext], allEdges[i], theFace, v );
+    //     if ( angle < M_PI ) // [-PI; PI]
+    //       convex += ( angle + M_PI ) / M_PI / M_PI;
+    //     isStraight[i] += cornerWgt * convex;
+    //   }
+    // }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Creates an EDGE from a sole branch of MA
+   */
+  //================================================================================
+
+  TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper&            theHelper,
+                              const SMESH_MAT2d::MedialAxis& theMA )
+  {
+    if ( theMA.getBranches().size() != 1 )
+      return TopoDS_Edge();
+
+    vector< gp_XY > uv;
+    theMA.getPoints( theMA.getBranches()[0], uv );
+    if ( uv.size() < 2 )
+      return TopoDS_Edge();
+
+    TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
+    Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
+
+    // cout << "from salome.geom import geomBuilder" << endl;
+    // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl;
+    Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size());
+    for ( size_t i = 0; i < uv.size(); ++i )
+    {
+      gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() );
+      points->SetValue( i+1, p );
+      //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl;
+    }
+
+    GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution());
+    interpol.Perform();
+    if ( !interpol.IsDone())
+      return TopoDS_Edge();
+
+    TopoDS_Edge branchEdge = BRepBuilderAPI_MakeEdge(interpol.Curve());
+    return branchEdge;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Returns a type of shape, to which a hypothesis used to mesh a given edge is assigned
+   */
+  //================================================================================
+
+  TopAbs_ShapeEnum getHypShape( SMESH_Mesh* mesh, const TopoDS_Shape& edge )
+  {
+    TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
+
+    SMESH_subMesh* sm = mesh->GetSubMesh( edge );
+    SMESH_Algo*  algo = sm->GetAlgo();
+    if ( !algo ) return shapeType;
+
+    const list <const SMESHDS_Hypothesis *> & hyps =
+      algo->GetUsedHypothesis( *mesh, edge, /*ignoreAuxiliary=*/true );
+    if ( hyps.empty() ) return shapeType;
+
+    TopoDS_Shape shapeOfHyp =
+      SMESH_MesherHelper::GetShapeOfHypothesis( hyps.front(), edge, mesh);
+
+    return SMESH_MesherHelper::GetGroupType( shapeOfHyp, /*woCompound=*/true);
+  }
+
+  //================================================================================
+  /*!
+   * \brief Discretize a sole branch of MA an returns parameters of divisions on MA
+   */
+  //================================================================================
+
+  bool divideMA( SMESH_MesherHelper&            theHelper,
+                 const SMESH_MAT2d::MedialAxis& theMA,
+                 const vector<TopoDS_Edge>&     theSinuEdges,
+                 const size_t                   theSinuSide0Size,
+                 SMESH_Algo*                    the1dAlgo,
+                 vector<double>&                theMAParams )
+  {
+    // check if all EDGEs of one size are meshed, then MA discretization is not needed
+    SMESH_Mesh* mesh = theHelper.GetMesh();
+    size_t nbComputedEdges[2] = { 0, 0 };
+    for ( size_t i = 1; i < theSinuEdges.size(); ++i )
+    {
+      bool isComputed = ( ! mesh->GetSubMesh( theSinuEdges[i] )->IsEmpty() );
+      nbComputedEdges[ i < theSinuSide0Size ] += isComputed;
+    }
+    if ( nbComputedEdges[0] == theSinuSide0Size ||
+         nbComputedEdges[1] == theSinuEdges.size() - theSinuSide0Size )
+      return true; // discretization is not needed
+
+
+    TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA );
+    if ( branchEdge.IsNull() )
+      return false;
+
+    // const char* file = "/misc/dn25/salome/eap/salome/misc/tmp/MAedge.brep";
+    // BRepTools::Write( branchEdge, file);
+    // cout << "Write " << file << endl;
+
+    // look for a most local hyps assigned to theSinuEdges
+    TopoDS_Edge edge = theSinuEdges[0];
+    int mostSimpleShape = (int) getHypShape( mesh, edge );
+    for ( size_t i = 1; i < theSinuEdges.size(); ++i )
+    {
+      int shapeType = (int) getHypShape( mesh, theSinuEdges[i] );
+      if ( shapeType > mostSimpleShape )
+        edge = theSinuEdges[i];
+    }
+
+    SMESH_Algo* algo = the1dAlgo;
+    if ( mostSimpleShape != TopAbs_SHAPE )
+    {
+      algo = mesh->GetSubMesh( edge )->GetAlgo();
+      SMESH_Hypothesis::Hypothesis_Status status;
+      if ( !algo->CheckHypothesis( *mesh, edge, status ))
+        algo = the1dAlgo;
+    }
+
+    TmpMesh tmpMesh;
+    tmpMesh.ShapeToMesh( branchEdge );
+    try {
+      mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes
+      if ( !algo->Compute( tmpMesh, branchEdge ))
+        return false;
+    }
+    catch (...) {
+      return false;
+    }
+    return SMESH_Algo::GetNodeParamOnEdge( tmpMesh.GetMeshDS(), branchEdge, theMAParams );
+  }
+
+  //================================================================================
+  /*!
+   * \brief Modifies division parameters on MA to make them coincide with projections
+   *        of VERTEXes to MA for a given pair of opposite EDGEs
+   *  \param [in] theEdgePairInd - index of the EDGE pair
+   *  \param [in] theDivPoints - the BranchPoint's dividing MA into parts each
+   *         corresponding to a unique pair of opposite EDGEs
+   *  \param [in,out] theMAParams - the MA division parameters to modify
+   *  \param [in,out] theParBeg - index of the 1st division point for the given EDGE pair
+   *  \param [in,out] theParEnd - index of the last division point for the given EDGE pair
+   *  \return bool - is OK
+   */
+  //================================================================================
+
+  bool getParamsForEdgePair( const size_t                              theEdgePairInd,
+                             const vector< SMESH_MAT2d::BranchPoint >& theDivPoints,
+                             const vector<double>&                     theMAParams,
+                             vector<double>&                           theSelectedMAParams)
+  {
+    if ( theDivPoints.empty() )
+    {
+      theSelectedMAParams = theMAParams;
+      return true;
+    }
+    if ( theEdgePairInd > theDivPoints.size() )
+      return false;
+
+    // TODO
+    return false;
+  }
+
+  //--------------------------------------------------------------------------------
+  // node or node parameter on EDGE
+  struct NodePoint
+  {
+    const SMDS_MeshNode* _node;
+    double               _u;
+    int                  _edgeInd; // index in theSinuEdges vector
+
+    NodePoint(const SMDS_MeshNode* n=0 ): _node(n), _u(0), _edgeInd(-1) {}
+    NodePoint(double u, size_t iEdge) : _node(0), _u(u), _edgeInd(iEdge) {}
+    NodePoint(const SMESH_MAT2d::BoundaryPoint& p) : _node(0), _u(p._param), _edgeInd(p._edgeIndex) {}
+  };
+
+  //================================================================================
+  /*!
+   * \brief Finds a VERTEX corresponding to a point on EDGE, which is also filled
+   *        with a node on the VERTEX, present or created
+   *  \param [in,out] theNodePnt - the node position on the EDGE
+   *  \param [in] theSinuEdges - the sinuous EDGEs
+   *  \param [in] theMeshDS - the mesh
+   *  \return bool - true if the \a theBndPnt is on VERTEX
+   */
+  //================================================================================
+
+  bool findVertex( NodePoint&                  theNodePnt,
+                   const vector<TopoDS_Edge>&  theSinuEdges,
+                   SMESHDS_Mesh*               theMeshDS)
+  {
+    if ( theNodePnt._edgeInd >= theSinuEdges.size() )
+      return false;
+
+    double f,l;
+    BRep_Tool::Range( theSinuEdges[ theNodePnt._edgeInd ], f,l );
+
+    TopoDS_Vertex V;
+    if      ( Abs( f - theNodePnt._u ))
+      V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
+    else if ( Abs( l - theNodePnt._u ))
+      V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false);
+
+    if ( !V.IsNull() )
+    {
+      theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS );
+      if ( !theNodePnt._node )
+      {
+        gp_Pnt p = BRep_Tool::Pnt( V );
+        theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() );
+        theMeshDS->SetNodeOnVertex( theNodePnt._node, V );
+      }
+      return true;
+    }
+    return false;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Add to the map of NodePoint's those on VERTEXes
+   *  \param [in,out] theHelper - the helper
+   *  \param [in] theMA - Medial Axis
+   *  \param [in] theDivPoints - projections of VERTEXes to MA
+   *  \param [in] theSinuEdges - the sinuous EDGEs
+   *  \param [in] theSideEdgeIDs - indices of sinuous EDGEs per side
+   *  \param [in] theIsEdgeComputed - is sinuous EGDE is meshed
+   *  \param [in,out] thePointsOnE - the map to fill
+   */
+  //================================================================================
+
+  bool projectVertices( SMESH_MesherHelper&                           theHelper,
+                        const SMESH_MAT2d::MedialAxis&                theMA,
+                        const vector< SMESH_MAT2d::BranchPoint >&     theDivPoints,
+                        const vector<TopoDS_Edge>&                    theSinuEdges,
+                        //const vector< int >                           theSideEdgeIDs[2],
+                        const vector< bool >&                         theIsEdgeComputed,
+                        map< double, pair< NodePoint, NodePoint > > & thePointsOnE)
+  {
+    if ( theDivPoints.empty() )
+      return true;
+
+    SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
+
+    double uMA;
+    SMESH_MAT2d::BoundaryPoint bp[2];
+    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+
+    for ( size_t i = 0; i < theDivPoints.size(); ++i )
+    {
+      if ( !branch.getParameter( theDivPoints[i], uMA ))
+        return false;
+      if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] ))
+        return false;
+
+      NodePoint np[2] = { NodePoint( bp[0] ),
+                          NodePoint( bp[1] ) };
+      bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ),
+                           findVertex( np[1], theSinuEdges, meshDS )};
+
+      map< double, pair< NodePoint, NodePoint > >::iterator u2NP =
+        thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first;
+
+      if ( isVertex[0] && isVertex[1] )
+        continue;
+
+      bool isOppComputed = theIsEdgeComputed[ np[ isVertex[0] ]._edgeInd ];
+
+      if ( !isOppComputed )
+        continue;
+
+      // a VERTEX is projected on a meshed EDGE; there are two options:
+      // - a projected point is joined with a closet node if a strip between this and neighbor
+      // projection is wide enough; joining is done by setting the same node to the BoundaryPoint
+      // - a neighbor projection is merged this this one if it too close; a node of deleted
+      // projection is set to the BoundaryPoint of this projection
+
+
+    }
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Divide the sinuous EDGEs by projecting the division point of Medial
+   *        Axis to the EGDEs
+   *  \param [in] theHelper - the helper
+   *  \param [in] theMA - the Medial Axis
+   *  \param [in] theMAParams - parameters of division points of \a theMA
+   *  \param [in] theSinuEdges - the EDGEs to make nodes on
+   *  \param [in] theSinuSide0Size - the number of EDGEs in the 1st sinuous side
+   *  \return bool - is OK or not
+   */
+  //================================================================================
+
+  bool computeSinuEdges( SMESH_MesherHelper&        theHelper,
+                         SMESH_MAT2d::MedialAxis&   theMA,
+                         vector<double>&            theMAParams,
+                         const vector<TopoDS_Edge>& theSinuEdges,
+                         const size_t               theSinuSide0Size)
+  {
+    if ( theMA.getBranches().size() != 1 )
+      return false;
+
+    // normalize theMAParams
+    for ( size_t i = 0; i < theMAParams.size(); ++i )
+      theMAParams[i] /= theMAParams.back();
+
+
+    SMESH_Mesh*     mesh = theHelper.GetMesh();
+    SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
+    double f,l;
+
+    vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() );
+    vector< int >                edgeIDs( theSinuEdges.size() );
+    vector< bool >            isComputed( theSinuEdges.size() );
+    //bool hasComputed = false;
+    for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+    {
+      curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l );
+      if ( !curves[i] )
+        return false;
+      SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
+      edgeIDs   [i] = sm->GetId();
+      isComputed[i] = ( !sm->IsEmpty() );
+      // if ( isComputed[i] )
+      //   hasComputed = true;
+    }
+
+    const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0];
+    SMESH_MAT2d::BoundaryPoint bp[2];
+
+    vector< std::size_t > edgeIDs1, edgeIDs2;
+    vector< SMESH_MAT2d::BranchPoint > divPoints;
+    branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints );
+    for ( size_t i = 0; i < edgeIDs1.size(); ++i )
+      if ( isComputed[ edgeIDs1[i]] &&
+           isComputed[ edgeIDs2[i]])
+        return false;
+
+    // map param on MA to parameters of nodes on a pair of theSinuEdges
+    typedef map< double, pair< NodePoint, NodePoint > > TMAPar2NPoints;
+    TMAPar2NPoints pointsOnE;
+    vector<double> maParams;
+
+    // compute params of nodes on EDGEs by projecting division points from MA
+    //const double tol = 1e-5 * theMAParams.back();
+    size_t iEdgePair = 0;
+    while ( iEdgePair < edgeIDs1.size() )
+    {
+      if ( isComputed[ edgeIDs1[ iEdgePair ]] ||
+           isComputed[ edgeIDs2[ iEdgePair ]])
+      {
+        // "projection" from one side to the other
+
+        size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0;
+        if ( !isComputed[ iEdgeComputed ])
+          ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair];
+
+        map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes
+        if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams ))
+          return false;
+
+        SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ];
+        SMESH_MAT2d::BranchPoint brp;
+        NodePoint npN, npB;
+        NodePoint& np0 = iSideComputed ? npB : npN;
+        NodePoint& np1 = iSideComputed ? npN : npB;
+
+        double maParam1st, maParamLast, maParam;
+        if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp ))
+            return false;
+        branch.getParameter( brp, maParam1st );
+        if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp ))
+            return false;
+        branch.getParameter( brp, maParamLast );
+
+        map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = --nodeParams.end();
+        TMAPar2NPoints::iterator pos, end = pointsOnE.end();
+        TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos;
+        for ( ++u2n; u2n != u2nEnd; ++u2n )
+        {
+          if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp ))
+            return false;
+          if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] ))
+            return false;
+          if ( !branch.getParameter( brp, maParam ))
+            return false;
+
+          npN = NodePoint( u2n->second );
+          npB = NodePoint( bndPnt );
+          pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 )));
+        }
+
+        // move iEdgePair forward
+        while ( iEdgePair < edgeIDs1.size() )
+          if ( edgeIDs1[ iEdgePair ] == bp[0]._edgeIndex &&
+               edgeIDs2[ iEdgePair ] == bp[1]._edgeIndex )
+            break;
+          else
+            ++iEdgePair;
+      }
+      else
+      {
+        // projection from MA
+        maParams.clear();
+        if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams ))
+          return false;
+
+        for ( size_t i = 1; i < maParams.size()-1; ++i )
+        {
+          if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] ))
+            return false;
+
+          pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]),
+                                                                                NodePoint(bp[1]))));
+        }
+      }
+      ++iEdgePair;
+    }
+
+    if ( !projectVertices( theHelper, theMA, divPoints, theSinuEdges, isComputed, pointsOnE ))
+      return false;
+
+    // create nodes
+    TMAPar2NPoints::iterator u2np = pointsOnE.begin();
+    for ( ; u2np != pointsOnE.end(); ++u2np )
+    {
+      NodePoint* np[2] = { & u2np->second.first, & u2np->second.second };
+      for ( int iSide = 0; iSide < 2; ++iSide )
+      {
+        if ( np[ iSide ]->_node ) continue;
+        size_t       iEdge = np[ iSide ]->_edgeInd;
+        double           u = np[ iSide ]->_u;
+        gp_Pnt           p = curves[ iEdge ]->Value( u );
+        np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() );
+        meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u );
+      }
+    }
+
+    // create mesh segments on EDGEs
+    theHelper.SetElementsOnShape( false );
+    TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() );
+    for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+    {
+      SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] );
+      if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 )
+        continue;
+
+      StdMeshers_FaceSide side( face, theSinuEdges[i], mesh,
+                                /*isFwd=*/true, /*skipMediumNodes=*/true );
+      vector<const SMDS_MeshNode*> nodes = side.GetOrderedNodes();
+      for ( size_t in = 1; in < nodes.size(); ++in )
+      {
+        const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false );
+        meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] );
+      }
+    }
+
+    // update sub-meshes on VERTEXes
+    for ( size_t i = 0; i < theSinuEdges.size(); ++i )
+    {
+      mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] ))
+        ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+      mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] ))
+        ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+    }
+
+    return true;
+  }
+
+  //================================================================================
+  /*!
+   * \brief Mesh short EDGEs
+   */
+  //================================================================================
+
+  bool computeShortEdges( SMESH_MesherHelper&        theHelper,
+                          const vector<TopoDS_Edge>& theShortEdges,
+                          SMESH_Algo*                the1dAlgo )
+  {
+    for ( size_t i = 0; i < theShortEdges.size(); ++i )
+    {
+      theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true );
+
+      SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] );
+      if ( sm->IsEmpty() )
+      {
+        try {
+          if ( !the1dAlgo->Compute( *theHelper.GetMesh(), theShortEdges[i] ))
+            return false;
+        }
+        catch (...) {
+          return false;
+        }
+        sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
+        if ( sm->IsEmpty() )
+          return false;
+      }
+    }
+    return true;
+  }
+
+  inline double area( const UVPtStruct& p1, const UVPtStruct& p2, const UVPtStruct& p3 )
+  {
+    gp_XY v1 = p2.UV() - p1.UV();
+    gp_XY v2 = p3.UV() - p1.UV();
+    return v2 ^ v1;
+  }
+
+  bool ellipticSmooth( FaceQuadStruct::Ptr quad, int nbLoops )
+  {
+    //nbLoops = 10;
+    if ( quad->uv_grid.empty() )
+      return true;
+
+    int nbhoriz  = quad->iSize;
+    int nbvertic = quad->jSize;
+
+    const double dksi = 0.5, deta = 0.5;
+    const double  dksi2 = dksi*dksi, deta2 = deta*deta;
+    double err = 0., g11, g22, g12;
+    int nbErr = 0;
+
+    FaceQuadStruct& q = *quad;
+    UVPtStruct pNew;
+
+    double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) );
+
+    for ( int iLoop = 0; iLoop < nbLoops; ++iLoop )
+    {
+      err = 0;
+      for ( int i = 1; i < nbhoriz - 1; i++ )
+        for ( int j = 1; j < nbvertic - 1; j++ )
+        {
+          g11 = ( (q.U(i,j+1) - q.U(i,j-1))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
+                  (q.V(i,j+1) - q.V(i,j-1))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/4;
+
+          g22 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i+1,j) - q.U(i-1,j))/dksi2 +
+                  (q.V(i+1,j) - q.V(i-1,j))*(q.V(i+1,j) - q.V(i-1,j))/deta2 )/4;
+
+          g12 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 +
+                  (q.V(i+1,j) - q.V(i-1,j))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/(4*dksi*deta);
+
+          pNew.u = dksi2/(2*(g11+g22)) * (g11*(q.U(i+1,j) + q.U(i-1,j))/dksi2 +
+                                          g22*(q.U(i,j+1) + q.U(i,j-1))/dksi2
+                                          - 0.5*g12*q.U(i+1,j+1) + 0.5*g12*q.U(i-1,j+1) +
+                                          - 0.5*g12*q.U(i-1,j-1) + 0.5*g12*q.U(i+1,j-1));
+
+          pNew.v = deta2/(2*(g11+g22)) * (g11*(q.V(i+1,j) + q.V(i-1,j))/deta2 +
+                                          g22*(q.V(i,j+1) + q.V(i,j-1))/deta2
+                                          - 0.5*g12*q.V(i+1,j+1) + 0.5*g12*q.V(i-1,j+1) +
+                                          - 0.5*g12*q.V(i-1,j-1) + 0.5*g12*q.V(i+1,j-1));
+
+          // if (( refArea * area( q.UVPt(i-1,j-1), q.UVPt(i,j-1), pNew ) > 0 ) &&
+          //     ( refArea * area( q.UVPt(i+1,j-1), q.UVPt(i+1,j), pNew ) > 0 ) &&
+          //     ( refArea * area( q.UVPt(i+1,j+1), q.UVPt(i,j+1), pNew ) > 0 ) &&
+          //     ( refArea * area( q.UVPt(i-1,j), q.UVPt(i-1,j-1), pNew ) > 0 ))
+          {
+            err += sqrt(( q.U(i,j) - pNew.u ) * ( q.U(i,j) - pNew.u ) +
+                        ( q.V(i,j) - pNew.v ) * ( q.V(i,j) - pNew.v ));
+            q.U(i,j) = pNew.u;
+            q.V(i,j) = pNew.v;
+          }
+          // else if ( ++nbErr < 10 )
+          // {
+          //   cout << i << ", " << j << endl;
+          //   cout << "x = ["
+          //        << "[ " << q.U(i-1,j-1) << ", " <<q.U(i,j-1) << ", " << q.U(i+1,j-1) << " ],"
+          //        << "[ " << q.U(i-1,j-0) << ", " <<q.U(i,j-0) << ", " << q.U(i+1,j-0) << " ],"
+          //        << "[ " << q.U(i-1,j+1) << ", " <<q.U(i,j+1) << ", " << q.U(i+1,j+1) << " ]]" << endl;
+          //   cout << "y = ["
+          //        << "[ " << q.V(i-1,j-1) << ", " <<q.V(i,j-1) << ", " << q.V(i+1,j-1) << " ],"
+          //        << "[ " << q.V(i-1,j-0) << ", " <<q.V(i,j-0) << ", " << q.V(i+1,j-0) << " ],"
+          //        << "[ " << q.V(i-1,j+1) << ", " <<q.V(i,j+1) << ", " << q.V(i+1,j+1) << " ]]" << endl<<endl;
+          // }
+        }
+
+      if ( err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) < 1e-6 )
+        break;
+    }
+    //cout << " ERR " << err / ( nbhoriz - 2 ) / ( nbvertic - 2 ) << endl;
+
+    return true;
+  }
+
+
+} // namespace
+
+//================================================================================
+/*!
+ * \brief Create quadrangle elements
+ *  \param [in] theHelper - the helper
+ *  \param [in] theFace - the face to mesh
+ *  \param [in] theSinuEdges - the sinuous EDGEs
+ *  \param [in] theShortEdges - the short EDGEs
+ *  \return bool - is OK or not
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::computeQuads( SMESH_MesherHelper&       theHelper,
+                                                       const TopoDS_Face&        theFace,
+                                                       const vector<TopoDS_Edge> theSinuEdges[2],
+                                                       const vector<TopoDS_Edge> theShortEdges[2])
+{
+  SMESH_Mesh* mesh = theHelper.GetMesh();
+  SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, theFace );
+  if ( !proxyMesh )
+    return false;
+
+  StdMeshers_Quadrangle_2D::myProxyMesh  = proxyMesh;
+  StdMeshers_Quadrangle_2D::myHelper     = &theHelper;
+  StdMeshers_Quadrangle_2D::myNeedSmooth = false;
+  StdMeshers_Quadrangle_2D::myCheckOri   = false;
+  StdMeshers_Quadrangle_2D::myQuadList.clear();
+
+  // fill FaceQuadStruct
+
+  list< TopoDS_Edge > side[4];
+  side[0].insert( side[0].end(), theShortEdges[0].begin(), theShortEdges[0].end() );
+  side[1].insert( side[1].end(), theSinuEdges[1].begin(),  theSinuEdges[1].end() );
+  side[2].insert( side[2].end(), theShortEdges[1].begin(), theShortEdges[1].end() );
+  side[3].insert( side[3].end(), theSinuEdges[0].begin(),  theSinuEdges[0].end() );
+
+  FaceQuadStruct::Ptr quad( new FaceQuadStruct );
+  quad->side.resize( 4 );
+  quad->face = theFace;
+  for ( int i = 0; i < 4; ++i )
+  {
+    quad->side[i] = StdMeshers_FaceSide::New( theFace, side[i], mesh, i < QUAD_TOP_SIDE,
+                                              /*skipMediumNodes=*/true, proxyMesh );
+  }
+  int nbNodesShort0 = quad->side[0].NbPoints();
+  int nbNodesShort1 = quad->side[2].NbPoints();
+
+  // compute UV of internal points
+  myQuadList.push_back( quad );
+  if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( quad ))
+    return false;
+
+  // elliptic smooth of internal points to get boundary cell normal to the boundary
+  ellipticSmooth( quad, 1 );
+
+  // create quadrangles
+  bool ok;
+  if ( nbNodesShort0 == nbNodesShort1 )
+    ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *mesh, theFace, quad );
+  else
+    ok = StdMeshers_Quadrangle_2D::computeTriangles( *mesh, theFace, quad );
+
+  StdMeshers_Quadrangle_2D::myProxyMesh.reset();
+  StdMeshers_Quadrangle_2D::myHelper = 0;
+  
+  return ok;
+}
+
+//================================================================================
+/*!
+ * \brief Generate quadrangle mesh
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
+                                                 const TopoDS_Shape& theShape)
+{
+  SMESH_MesherHelper helper( theMesh );
+  helper.SetSubShape( theShape );
+
+  TopoDS_Face F = TopoDS::Face( theShape );
+  if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
+
+  list< TopoDS_Edge > edges;
+  list< int > nbEdgesInWire;
+  int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
+
+  vector< TopoDS_Edge > sinuSide[2], shortSide[2];
+  if ( nbWire == 1 && getSinuousEdges( helper, F, edges, sinuSide, shortSide ))
+  {
+    vector< TopoDS_Edge > sinuEdges  = sinuSide[0];
+    sinuEdges.insert( sinuEdges.end(), sinuSide[1].begin(), sinuSide[1].end() );
+    if ( sinuEdges.size() > 2 )
+      return error(COMPERR_BAD_SHAPE, "Not yet supported case" );
+
+    double minSegLen = getMinSegLen( helper, sinuEdges );
+    SMESH_MAT2d::MedialAxis ma( F, sinuEdges, minSegLen, /*ignoreCorners=*/true );
+
+    if ( !_regular1D )
+      _regular1D = new Algo1D( _studyId, _gen );
+    _regular1D->SetSegmentLength( minSegLen );
+
+    vector<double> maParams;
+    if ( ! divideMA( helper, ma, sinuEdges, sinuSide[0].size(), _regular1D, maParams ))
+      return error(COMPERR_BAD_SHAPE);
+
+    if ( !computeShortEdges( helper, shortSide[0], _regular1D ) ||
+         !computeShortEdges( helper, shortSide[1], _regular1D ))
+      return error("Failed to mesh short edges");
+
+    if ( !computeSinuEdges( helper, ma, maParams, sinuEdges, sinuSide[0].size() ))
+      return error("Failed to mesh sinuous edges");
+
+    return computeQuads( helper, F, sinuSide, shortSide );
+  }
+
+  return error(COMPERR_BAD_SHAPE, "Not implemented so far");
+}
+
+//================================================================================
+/*!
+ * \brief Predict nb of elements
+ */
+//================================================================================
+
+bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh &         theMesh,
+                                                  const TopoDS_Shape & theShape,
+                                                  MapShapeNbElems&     theResMap)
+{
+  return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap);
+}
+
diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx
new file mode 100644
index 000000000..df8a84a4c
--- /dev/null
+++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx
@@ -0,0 +1,65 @@
+// Copyright (C) 2007-2015  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
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : StdMeshers_QuadFromMedialAxis_1D2D.hxx
+// Created   : Wed Jun  3 17:22:35 2015
+// Author    : Edward AGAPOV (eap)
+
+
+#ifndef __StdMeshers_QuadFromMedialAxis_1D2D_HXX__
+#define __StdMeshers_QuadFromMedialAxis_1D2D_HXX__
+
+#include "StdMeshers_Quadrangle_2D.hxx"
+
+#include <vector>
+
+/*!
+ * \brief Quadrangle mesher using Medial Axis
+ */
+class STDMESHERS_EXPORT StdMeshers_QuadFromMedialAxis_1D2D: public StdMeshers_Quadrangle_2D
+{
+ public:
+  StdMeshers_QuadFromMedialAxis_1D2D(int hypId, int studyId, SMESH_Gen* gen);
+  virtual ~StdMeshers_QuadFromMedialAxis_1D2D();
+
+  virtual bool CheckHypothesis(SMESH_Mesh&         aMesh,
+                               const TopoDS_Shape& aShape,
+                               Hypothesis_Status&  aStatus);
+
+  virtual bool Compute(SMESH_Mesh&         aMesh,
+                       const TopoDS_Shape& aShape);
+
+  virtual bool Evaluate(SMESH_Mesh &         aMesh,
+                        const TopoDS_Shape & aShape,
+                        MapShapeNbElems&     aResMap);
+
+ private:
+
+  bool computeQuads( SMESH_MesherHelper&            theHelper,
+                     const TopoDS_Face&             theFace,
+                     const std::vector<TopoDS_Edge> theSinuEdges[2],
+                     const std::vector<TopoDS_Edge> theShortEdges[2]);
+
+  class Algo1D;
+  Algo1D* _regular1D;
+};
+
+#endif
diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
index af3d51d76..0e35c560e 100644
--- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
+++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx
@@ -1599,7 +1599,7 @@ void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int nu
 
 //================================================================================
 /*!
- * \brief Rotate sides of a quad by given nb of quartes
+ * \brief Rotate sides of a quad CCW by given nb of quartes
  *  \param nb  - number of rotation quartes
  *  \param ori - to keep orientation of sides as in an unit quad or not
  *  \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to
@@ -1611,6 +1611,8 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid )
 {
   if ( nb == 0 ) return;
 
+  nb = nb % NB_QUAD_SIDES;
+
   vector< Side > newSides( side.size() );
   vector< Side* > sidePtrs( side.size() );
   for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i)
@@ -1640,7 +1642,33 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid )
   }
   newSides.swap( side );
 
-  uv_grid.clear();
+  if ( keepGrid && !uv_grid.empty() )
+  {
+    if ( nb == 2 ) // "PI"
+    {
+      std::reverse( uv_grid.begin(), uv_grid.end() );
+    }
+    else
+    {
+      FaceQuadStruct newQuad;
+      newQuad.uv_grid.resize( uv_grid.size() );
+      newQuad.iSize = jSize;
+      newQuad.jSize = iSize;
+      int i, j, iRev, jRev;
+      int *iNew = ( nb == 1 ) ? &jRev : &j;
+      int *jNew = ( nb == 1 ) ? &i : &iRev;
+      for ( i = 0, iRev = iSize-1; i < iSize; ++i, --iRev )
+        for ( j = 0, jRev = jSize-1; j < jSize; ++j, --jRev )
+          newQuad.UVPt( *iNew, *jNew ) = UVPt( i, j );
+
+      std::swap( iSize, jSize );
+      std::swap( uv_grid, newQuad.uv_grid );
+    }
+  }
+  else
+  {
+    uv_grid.clear();
+  }
 }
 
 //=======================================================================
diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx
index 8687bbfda..6a06c9a46 100644
--- a/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx
+++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx
@@ -119,6 +119,8 @@ struct FaceQuadStruct
 
   FaceQuadStruct ( const TopoDS_Face& F = TopoDS_Face(), const std::string& nm="main" );
   UVPtStruct& UVPt( int i, int j ) { return uv_grid[ i + j * iSize ]; }
+  double&        U( int i, int j ) { return UVPt( i, j ).u; }
+  double&        V( int i, int j ) { return UVPt( i, j ).v; }
   void  shift    ( size_t nb, bool keepUnitOri, bool keepGrid=false );
   int & nbNodeOut( int iSide ) { return side[ iSide ].nbNodeOut; }
   bool  findCell ( const gp_XY& uv, int & i, int & j );
diff --git a/src/StdMeshersGUI/StdMeshers_images.ts b/src/StdMeshersGUI/StdMeshers_images.ts
index b94521ab3..e744cf34e 100644
--- a/src/StdMeshersGUI/StdMeshers_images.ts
+++ b/src/StdMeshersGUI/StdMeshers_images.ts
@@ -143,6 +143,10 @@
             <source>ICON_SMESH_TREE_ALGO_RadialQuadrangle_1D2D</source>
             <translation>mesh_tree_algo_radial_quadrangle_1D2D.png</translation>
         </message>
+	<message>
+            <source>ICON_SMESH_TREE_ALGO_QuadFromMedialAxis_1D2D</source>
+            <translation>mesh_tree_algo_quad.png</translation>
+        </message>
 	<message>
             <source>ICON_SMESH_TREE_ALGO_Prism_3D</source>
             <translation>mesh_tree_algo_prism.png</translation>
diff --git a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx
index 2af7cca75..cf497681c 100644
--- a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx
+++ b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx
@@ -32,6 +32,8 @@
 #include "Utils_CorbaException.hxx"
 #include "utilities.h"
 
+#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx"
+
 using namespace std;
 
 //=============================================================================
@@ -98,3 +100,39 @@ CORBA::Boolean StdMeshers_Quadrangle_2D_i::IsApplicable( const TopoDS_Shape &S,
   return ::StdMeshers_Quadrangle_2D::IsApplicable( S, toCheckAll );
 }
 
+//=============================================================================
+/*!
+ *  StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i
+ *
+ *  Constructor
+ */
+//=============================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D_i::
+StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA,
+                                      int                     theStudyId,
+                                      ::SMESH_Gen*            theGenImpl )
+  : SALOME::GenericObj_i( thePOA ),
+    SMESH_Hypothesis_i( thePOA ),
+    SMESH_Algo_i( thePOA ),
+    SMESH_2D_Algo_i( thePOA )
+{
+  MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i" );
+  myBaseImpl = new ::StdMeshers_QuadFromMedialAxis_1D2D( theGenImpl->GetANewId(),
+                                                         theStudyId,
+                                                         theGenImpl );
+}
+
+//=============================================================================
+/*!
+ *  StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i
+ *
+ *  Destructor
+ *
+ */
+//=============================================================================
+
+StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i()
+{
+  MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i" );
+}
diff --git a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx
index 4fa55b135..145d31e3c 100644
--- a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx
+++ b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx
@@ -62,4 +62,27 @@ class STDMESHERS_I_EXPORT StdMeshers_Quadrangle_2D_i:
   static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll);
 };
 
+// ======================================================
+// Quadrangle (Medial Axis Projection) 2d algorithm
+// ======================================================
+class STDMESHERS_I_EXPORT StdMeshers_QuadFromMedialAxis_1D2D_i:
+  public virtual POA_StdMeshers::StdMeshers_QuadFromMedialAxis_1D2D,
+  public virtual SMESH_2D_Algo_i
+{
+ public:
+  // Constructor
+  StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA,
+                                        int                     theStudyId,
+                                        ::SMESH_Gen*            theGenImpl );
+
+  // Destructor
+  virtual ~StdMeshers_QuadFromMedialAxis_1D2D_i();
+
+  // Get implementation
+  //::StdMeshers_Quadrangle_2D* GetImpl();
+
+  // Return true if the algorithm is applicable to a shape
+  //static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll);
+};
+
 #endif
diff --git a/src/StdMeshers_I/StdMeshers_i.cxx b/src/StdMeshers_I/StdMeshers_i.cxx
index c7ae1a76a..3e9eb70b7 100644
--- a/src/StdMeshers_I/StdMeshers_i.cxx
+++ b/src/StdMeshers_I/StdMeshers_i.cxx
@@ -216,6 +216,8 @@ STDMESHERS_I_EXPORT
 #endif
     else if (strcmp(aHypName, "Quadrangle_2D") == 0)
       aCreator = new StdHypothesisCreator_i<StdMeshers_Quadrangle_2D_i, StdMeshers_Quadrangle_2D_i>;
+    else if (strcmp(aHypName, "QuadFromMedialAxis_1D2D") == 0)
+      aCreator = new StdHypothesisCreator_i<StdMeshers_QuadFromMedialAxis_1D2D_i, StdMeshers_Quadrangle_2D_i>;
     else if (strcmp(aHypName, "Hexa_3D") == 0)
       aCreator = new StdHypothesisCreator_i<StdMeshers_Hexa_3D_i, StdMeshers_Hexa_3D_i>;
     else if (strcmp(aHypName, "Projection_1D") == 0)