// Copyright (C) 2007-2013  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.
//
// 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_Adaptive1D.cxx
//  Module : SMESH
//
#include "StdMeshers_Adaptive1D.hxx"

#include "SMESH_Gen.hxx"
#include "SMESH_Mesh.hxx"
#include "SMESH_MesherHelper.hxx"
#include "SMESH_Octree.hxx"
#include "SMESH_subMesh.hxx"
#include "SMESH_HypoFilter.hxx"

#include <Utils_SALOME_Exception.hxx>

#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepMesh_IncrementalMesh.hxx>
#include <BRep_Tool.hxx>
#include <Bnd_B3d.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <GeomAdaptor_Curve.hxx>
#include <Geom_Curve.hxx>
#include <Poly_Array1OfTriangle.hxx>
#include <Poly_Triangulation.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopLoc_Location.hxx>
#include <TopTools_IndexedMapOfShape.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Vertex.hxx>
#include <gp_Lin.hxx>
#include <gp_Pnt.hxx>

#include <limits>
#include <vector>
#include <set>

using namespace std;

namespace // internal utils
{
  //================================================================================
  /*!
   * \brief Working data of an EDGE
   */
  struct EdgeData
  {
    struct ProbePnt
    {
      gp_Pnt myP;
      double myU;
      double mySegSize;
      ProbePnt( gp_Pnt p, double u, double sz=1e100 ): myP( p ), myU( u ), mySegSize( sz ) {}
    };
    BRepAdaptor_Curve myC3d;
    double            myLength;
    list< ProbePnt >  myPoints;

    typedef list< ProbePnt >::iterator TPntIter;
    void AddPoint( TPntIter where, double u )
    {
      myPoints.insert( where, ProbePnt( myC3d.Value( u ), u ));
    }
    const ProbePnt& First() const { return myPoints.front(); }
    const ProbePnt& Last()  const { return myPoints.back(); }
    const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
  };

  //================================================================================
  /*!
   * \brief Bnd_B3d with access to its center and half-size
   */
  struct BBox : public Bnd_B3d
  {
    gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); }
    gp_XYZ HSize()  const { return gp_XYZ( myHSize[0],  myHSize[1],  myHSize[2]  ); }
    double Size()   const { return 2 * myHSize[0]; }
  };
  //================================================================================
  /*!
   * \brief Octree of local segment size
   */
  class SegSizeTree : public SMESH_Octree
  {
    double mySegSize; // segment size

    // structure holding some common parameters of SegSizeTree
    struct _CommonData : public SMESH_TreeLimit
    {
      double myGrading, myMinSize, myMaxSize;
    };
    _CommonData* getData() const { return (_CommonData*) myLimit; }

    SegSizeTree(double size): SMESH_Octree(), mySegSize(size)
    {
      allocateChildren();
    }
    void allocateChildren()
    {
      myChildren = new SMESH_Octree::TBaseTree*[nbChildren()];
      for ( int i = 0; i < nbChildren(); ++i )
        myChildren[i] = NULL;
    }
    virtual box_type* buildRootBox() { return 0; }
    virtual SegSizeTree* newChild() const { return 0; }
    virtual void buildChildrenData() {}

  public:

    SegSizeTree( Bnd_B3d & bb, double grading, double mixSize, double maxSize);
    void   SetSize( const gp_Pnt& p, double size );
    double SetSize( const gp_Pnt& p1, const gp_Pnt& p2 );
    double GetSize( const gp_Pnt& p ) const;
    const BBox* GetBox() const { return (BBox*) getBox(); }
    double GetMinSize() { return getData()->myMinSize; }
  };
  //================================================================================
  /*!
   * \brief Data of triangle used to locate it in an octree and to find distance
   *        to a point
   */
  struct Triangle
  {
    Bnd_B3d myBox;
    // data for DistToProjection()
    gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
    double myInvDet, myMaxSize2;

    void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 ); 
    bool DistToProjection( const gp_Pnt& p, double& dist ) const;
  };
  //================================================================================
  /*!
   * \brief Element data held by ElementBndBoxTree + algorithm computing a distance
   *        from a point to element
   */
  class ElementBndBoxTree;
  struct ElemTreeData : public SMESH_TreeLimit
  {
    virtual const Bnd_B3d* GetBox(int elemID) const = 0;
  };
  struct TriaTreeData : public ElemTreeData
  {
    vector< Triangle >           myTrias;
    double                       myFaceTol;
    double                       myTriasDeflection;
    const ElementBndBoxTree*     myTree;
    const Poly_Array1OfTriangle* myPolyTrias;
    const TColgp_Array1OfPnt*    myNodes;
    bool                         myOwnNodes;

    TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
    ~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; }
    virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
    void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const;
    double GetMinDistInSphere(const gp_Pnt& p,
                              const double  radius,
                              const bool    projectedOnly,
                              const gp_Pnt* avoidP=0) const;
  };
  //================================================================================
  /*!
   * \brief Octree of triangles or segments
   */
  class ElementBndBoxTree : public SMESH_Octree
  {
  public:
    ElementBndBoxTree(const TopoDS_Face& face);
    void GetElementsInSphere( const gp_XYZ& center,
                              const double  radius, std::set<int> & foundElemIDs) const;
    ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
    TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }

  protected:
    ElementBndBoxTree() {}
    SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
    void          buildChildrenData();
    Bnd_B3d*      buildRootBox();
  private:
    vector< int > _elementIDs;
  };
  //================================================================================
  /*!
   * \brief BRepMesh_IncrementalMesh with access to its protected Bnd_Box
   */
  struct IncrementalMesh : public BRepMesh_IncrementalMesh
  {
    IncrementalMesh(const TopoDS_Shape& shape,
                    const Standard_Real deflection,
                    const bool          relative):
      BRepMesh_IncrementalMesh( shape, deflection, relative )
    {
    }
    Bnd_B3d GetBox() const
    {
      Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
      myBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
      Bnd_B3d bb;
      bb.Add( gp_XYZ( TXmin, TYmin, TZmin ));
      bb.Add( gp_XYZ( TXmax, TYmax, TZmax ));
      return bb;
    }
  };

  //================================================================================
  /*!
   * \brief Initialize TriaTreeData
   */
  //================================================================================

  TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree )
    : myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false), myTriasDeflection(0)
  {
    TopLoc_Location loc;
    Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
    if ( !tr.IsNull() )
    {
      myFaceTol         = SMESH_MesherHelper::MaxTolerance( face );
      myTree            = triaTree;
      myNodes           = & tr->Nodes();
      myPolyTrias       = & tr->Triangles();
      myTriasDeflection = tr->Deflection();
      if ( !loc.IsIdentity() ) // transform nodes if necessary
      {
        TColgp_Array1OfPnt* trsfNodes = new TColgp_Array1OfPnt( myNodes->Lower(), myNodes->Upper() );
        trsfNodes->Assign( *myNodes );
        myNodes    = trsfNodes;
        myOwnNodes = true;
        const gp_Trsf& trsf = loc;
        for ( int i = trsfNodes->Lower(); i <= trsfNodes->Upper(); ++i )
          trsfNodes->ChangeValue(i).Transform( trsf );
      }
      myTrias.resize( myPolyTrias->Length() );
      Standard_Integer n1,n2,n3;
      for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
      {
        myPolyTrias->Value( i ).Get( n1,n2,n3 );
        myTrias[ i-1 ].Init( myNodes->Value( n1 ),
                             myNodes->Value( n2 ),
                             myNodes->Value( n3 ));
      }
      // TODO: mark triangles with nodes on VERTEXes to
      // less frequently compare with avoidPnt in GetMinDistInSphere()
      //
      // Handle(Poly_PolygonOnTriangulation) polygon =
      //   BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
      // if ( polygon.IsNull() /*|| !pologon.HasParameters()*/ )
      //   continue;
      // Handle(TColStd_Array1OfInteger) nodeIDs = polygon->Nodes();
    }
  }
  //================================================================================
  /*!
   * \brief Set size of segments by size of triangles
   */
  //================================================================================

  void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const
  {
    if ( myTriasDeflection <= std::numeric_limits<double>::min() )
      return;
    const double factor = deflection / myTriasDeflection;

    Standard_Integer n1,n2,n3;
    gp_Pnt p[4];
    double a[3];
    for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
    {
      // compute minimal altitude of a triangle
      myPolyTrias->Value( i ).Get( n1,n2,n3 );
      p[0] = myNodes->Value( n1 );
      p[1] = myNodes->Value( n2 );
      p[2] = myNodes->Value( n3 );
      p[3] = p[0];
      a[0] = p[0].Distance( p[1] );
      a[1] = p[1].Distance( p[2] );
      a[2] = p[2].Distance( p[3] );
      double maxSide = Max( a[0], Max( a[1], a[2] ));
      double s       = 0.5 * ( a[0] + a[1] + a[2] );
      double area    = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
      double sz      = 2 * area / maxSide; // minimal altitude
      for ( int j = 0; j < 3; ++j )
      {
        int nb = 2 * int( a[j] / sz + 0.5 );
        for ( int k = 1; k <= nb; ++k )
        {
          double r = double( k ) / nb;
          sizeTree.SetSize( r * p[j].XYZ() + ( 1-r ) * p[j+1].XYZ(), sz * factor );
        }
      }
      sizeTree.SetSize(( p[0].XYZ() + p[1].XYZ() + p[2].XYZ()) / 3., sz * factor );
      //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl;
    }
  }
  //================================================================================
  /*!
   * \brief Return minimal distance from a given point to a trinangle but not more
   *        distant than a given radius. Triangles with a node at avoidPnt are ignored.
   *        If projectedOnly, 
   */
  //================================================================================

  double TriaTreeData::GetMinDistInSphere(const gp_Pnt& p,
                                          const double  radius,
                                          const bool    projectedOnly,
                                          const gp_Pnt* avoidPnt) const
  {
    double minDist2 = 1e100;
    const double tol2 = myFaceTol * myFaceTol;

    std::set<int> foundElemIDs;
    myTree->GetElementsInSphere( p.XYZ(), radius, foundElemIDs );

    std::set<int>::iterator id = foundElemIDs.begin();
    Standard_Integer n[ 3 ];
    for ( ; id != foundElemIDs.end(); ++id )
    {
      double d, minD2 = minDist2;
      bool avoidTria = false;
      myPolyTrias->Value( *id+1 ).Get( n[0],n[1],n[2] );
      for ( int i = 0; i < 3; ++i )
      {
        const gp_Pnt& pn = myNodes->Value(n[i]);
        if ( avoidTria = ( avoidPnt && pn.SquareDistance(*avoidPnt) <= tol2 ))
          break;
        if ( !projectedOnly )
          minD2 = Min( minD2, pn.SquareDistance( p ));
      }
      if ( !avoidTria )
      {
        const Triangle& t = myTrias[ *id ];
        if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
          minD2 = Min( minD2, d*d );
        minDist2 = Min( minDist2, minD2 );
      }
    }
    return sqrt( minDist2 );
  }
  //================================================================================
  /*!
   * \brief Prepare Triangle data
   */
  //================================================================================

  void Triangle::Init( const gp_Pnt& p1, const gp_Pnt& p2, const gp_Pnt& p3 )
  {
    myBox.Add( p1 );
    myBox.Add( p2 );
    myBox.Add( p3 );
    myN0 = p1.XYZ();
    myEdge1 = p2.XYZ() - myN0;
    myEdge2 = p3.XYZ() - myN0;
    myNorm = myEdge1 ^ myEdge2;
    double normSize = myNorm.Modulus();
    if ( normSize > std::numeric_limits<double>::min() )
    {
      myNorm /= normSize;
      myPVec = myNorm ^ myEdge2;
      myInvDet = 1. / ( myEdge1 * myPVec );
    }
    else
    {
      myInvDet = 0.;
    }
    myMaxSize2 = Max( p2.SquareDistance( p3 ),
                      Max( myEdge2.SquareModulus(), myEdge1.SquareModulus() ));
  }
  //================================================================================
  /*!
   * \brief Compute distance from a point to the triangle. Return false if the point
   *        is not projected inside the triangle
   */
  //================================================================================

  bool Triangle::DistToProjection( const gp_Pnt& p, double& dist ) const
  {
    if ( myInvDet == 0 )
      return false; // degenerated triangle

    /* distance from n0 to the point */
    gp_XYZ tvec = p.XYZ() - myN0;

    /* calculate U parameter and test bounds */
    double u = ( tvec * myPVec ) * myInvDet;
    if (u < 0.0 || u > 1.0)
      return false; // projected outside the triangle

    /* calculate V parameter and test bounds */
    gp_XYZ qvec = tvec ^ myEdge1;
    double v = ( myNorm * qvec) * myInvDet;
    if ( v < 0.0 || u + v > 1.0 )
      return false; // projected outside the triangle

    dist = ( myEdge2 * qvec ) * myInvDet;
    return true;
  }

  //================================================================================
  /*!
   * \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE
   */
  //================================================================================

  ElementBndBoxTree::ElementBndBoxTree(const TopoDS_Face& face)
    :SMESH_Octree()
  {
    TriaTreeData* data = new TriaTreeData( face, this );
    data->myMaxLevel = 5;
    myLimit = data;

    if ( !data->myTrias.empty() )
    {
      for ( size_t i = 0; i < data->myTrias.size(); ++i )
        _elementIDs.push_back( i );

      compute();
    }
  }
  //================================================================================
  /*!
   * \brief Return the maximal box
   */
  //================================================================================

  Bnd_B3d* ElementBndBoxTree::buildRootBox()
  {
    Bnd_B3d* box = new Bnd_B3d;
    ElemTreeData* data = GetElemData();
    for ( int i = 0; i < _elementIDs.size(); ++i )
      box->Add( *data->GetBox( _elementIDs[i] ));
    return box;
  }

  //================================================================================
  /*!
   * \brief Redistrubute element boxes among children
   */
  //================================================================================

  void ElementBndBoxTree::buildChildrenData()
  {
    ElemTreeData* data = GetElemData();
    for ( int i = 0; i < _elementIDs.size(); ++i )
    {
      const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
      for (int j = 0; j < 8; j++)
        if ( !elemBox->IsOut( *myChildren[j]->getBox() ))
          ((ElementBndBoxTree*)myChildren[j])->_elementIDs.push_back( _elementIDs[i] );
    }
    SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory

    const int theMaxNbElemsInLeaf = 7;

    for (int j = 0; j < 8; j++)
    {
      ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
      if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
        child->myIsLeaf = true;
    }
  }
  //================================================================================
  /*!
   * \brief Return elements from leaves intersecting the sphere
   */
  //================================================================================

  void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ&   center,
                                               const double    radius,
                                               std::set<int> & foundElemIDs) const
  {
    if ( const box_type* box = getBox() )
    {
      if ( box->IsOut( center, radius ))
        return;

      if ( isLeaf() )
      {
        ElemTreeData* data = GetElemData();
        for ( int i = 0; i < _elementIDs.size(); ++i )
          if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
            foundElemIDs.insert( _elementIDs[i] );
      }
      else
      {
        for (int i = 0; i < 8; i++)
          ((ElementBndBoxTree*) myChildren[i])->GetElementsInSphere( center, radius, foundElemIDs );
      }
    }
  }

  //================================================================================
  /*!
   * \brief Constructor of SegSizeTree
   *  \param [in,out] bb - bounding box enclosing all EDGEs to discretize 
   *  \param [in] grading - factor to get max size of the neighbour segment by 
   *         size of a current one.
   */
  //================================================================================

  SegSizeTree::SegSizeTree( Bnd_B3d & bb, double grading, double minSize, double maxSize )
    : SMESH_Octree( new _CommonData() )
  {
    // make cube myBox from the box bb
    gp_XYZ pmin = bb.CornerMin(), pmax = bb.CornerMax();
    double maxBoxHSize = 0.5 * Max( pmax.X()-pmin.X(), Max( pmax.Y()-pmin.Y(), pmax.Z()-pmin.Z() ));
    maxBoxHSize *= 1.01;
    bb.SetHSize( gp_XYZ( maxBoxHSize, maxBoxHSize, maxBoxHSize ));
    myBox = new box_type( bb );

    mySegSize = Min( 2 * maxBoxHSize, maxSize );

    getData()->myGrading = grading;
    getData()->myMinSize = Max( minSize, 2*maxBoxHSize / 1.e6 );
    getData()->myMaxSize = maxSize;
    allocateChildren();
  }

  //================================================================================
  /*!
   * \brief Set segment size at a given point
   */
  //================================================================================

  void SegSizeTree::SetSize( const gp_Pnt& p, double size )
  {
    // check if the point is out of the largest cube
    SegSizeTree* root = this;
    while ( root->myFather )
      root = (SegSizeTree*) root->myFather;
    if ( root->getBox()->IsOut( p.XYZ() ))
      return;

    // keep size whthin the valid range
    size = Max( size, getData()->myMinSize );
    //size = Min( size, getData()->myMaxSize );

    // find an existing leaf at the point
    SegSizeTree* leaf = (SegSizeTree*) root;
    int iChild;
    while ( true )
    {
      iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
      if ( leaf->myChildren[ iChild ] )
        leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
      else
        break;
    }
    // don't increase the current size
    if ( leaf->mySegSize <= 1.1 * size )
      return;

    // split the found leaf until its box size is less than the given size
    const double rootSize = root->GetBox()->Size();
    while ( leaf->GetBox()->Size() > size )
    {
      const BBox* bb = leaf->GetBox();
      iChild   = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), bb->Center() );
      SegSizeTree* newLeaf = new SegSizeTree( bb->Size() / 2 );
      leaf->myChildren[iChild] = newLeaf;
      newLeaf->myFather = leaf;
      newLeaf->myLimit  = leaf->myLimit;
      newLeaf->myLevel  = leaf->myLevel + 1;
      newLeaf->myBox    = leaf->newChildBox( iChild );
      newLeaf->myBox->Enlarge( rootSize * 1e-10 );
      //newLeaf->myIsLeaf = ( newLeaf->mySegSize <= size );
      leaf = newLeaf;
    }
    leaf->mySegSize = size;

    // propagate increased size out from the leaf
    double boxSize = leaf->GetBox()->Size();
    double sizeInc = size + boxSize * getData()->myGrading;
    for ( int iDir = 1; iDir <= 3; ++iDir )
    {
      gp_Pnt outPnt = p;
      outPnt.SetCoord( iDir, p.Coord( iDir ) + boxSize );
      SetSize( outPnt, sizeInc );
      outPnt.SetCoord( iDir, p.Coord( iDir ) - boxSize );
      SetSize( outPnt, sizeInc );
    }
  }
  //================================================================================
  /*!
   * \brief Set size of a segment given by two end points
   */
  //================================================================================

  double SegSizeTree::SetSize( const gp_Pnt& p1, const gp_Pnt& p2 )
  {
    const double size = p1.Distance( p2 );
    gp_XYZ p = 0.5 * ( p1.XYZ() + p2.XYZ() );
    SetSize( p, size );
    SetSize( p1, size );
    SetSize( p2, size );
    //cout << "SetSize " << p1.Distance( p2 ) << " at " << p.X() <<", "<< p.Y()<<", "<<p.Z()<< endl;
    return GetSize( p );
  }

  //================================================================================
  /*!
   * \brief Return segment size at a point
   */
  //================================================================================

  double SegSizeTree::GetSize( const gp_Pnt& p ) const
  {
    const SegSizeTree* leaf = this;
    while ( true )
    {
      int iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
      if ( leaf->myChildren[ iChild ] )
        leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
      else
        return leaf->mySegSize;
    }
    return mySegSize; // just to return anything
  }

  //================================================================================
  /*!
   * \brief Evaluate curve deflection between two points
   * \param theCurve - the curve
   * \param theU1 - the parameter of the first point
   * \param theU2 - the parameter of the second point
   * \retval double - square deflection value
   */
  //================================================================================

  double deflection2(const BRepAdaptor_Curve & theCurve,
                     double                    theU1,
                     double                    theU2)
  {
    // line between theU1 and theU2
    gp_Pnt p1 = theCurve.Value( theU1 ), p2 = theCurve.Value( theU2 );
    gp_Lin segment( p1, gp_Vec( p1, p2 ));

    // evaluate square distance of theCurve from the segment
    Standard_Real dist2 = 0;
    const int nbPnt = 5;
    const double step = ( theU2 - theU1 ) / nbPnt;
    while (( theU1 += step ) < theU2 )
      dist2 = Max( dist2, segment.SquareDistance( theCurve.Value( theU1 )));

    return dist2;
  }

} // namespace

//=======================================================================
//function : StdMeshers_Adaptive1D
//purpose  : Constructor
StdMeshers_Adaptive1D::StdMeshers_Adaptive1D(int         hypId,
                                             int         studyId,
                                             SMESH_Gen * gen)
  :SMESH_Hypothesis(hypId, studyId, gen)
{
  myMinSize       = 1e-10;
  myMaxSize       = 1e+10;
  myDeflection    = 1e-2;
  myAlgo          = NULL;
  _name           = "Adaptive1D";
  _param_algo_dim = 1; // is used by SMESH_Regular_1D
}
//=======================================================================
//function : ~StdMeshers_Adaptive1D
//purpose  : Destructor
StdMeshers_Adaptive1D::~StdMeshers_Adaptive1D()
{
  delete myAlgo; myAlgo = NULL;
}
//=======================================================================
//function : SetDeflection
//purpose  : 
void StdMeshers_Adaptive1D::SetDeflection(double value)
  throw(SALOME_Exception)
{
  if (value <= std::numeric_limits<double>::min() )
    throw SALOME_Exception("Deflection must be greater that zero");
  if (myDeflection != value)
  {
    myDeflection = value;
    NotifySubMeshesHypothesisModification();
  }
}
//=======================================================================
//function : SetMinSize
//purpose  : Sets minimal allowed segment length
void StdMeshers_Adaptive1D::SetMinSize(double minSize)
  throw(SALOME_Exception)
{
  if (minSize <= std::numeric_limits<double>::min() )
    throw SALOME_Exception("Min size must be greater that zero");

  if (myMinSize != minSize )
  {
    myMinSize = minSize;
    NotifySubMeshesHypothesisModification();
  }
}
//=======================================================================
//function : SetMaxSize
//purpose  : Sets maximal allowed segment length
void StdMeshers_Adaptive1D::SetMaxSize(double maxSize)
  throw(SALOME_Exception)
{
  if (maxSize <= std::numeric_limits<double>::min() )
    throw SALOME_Exception("Max size must be greater that zero");

  if (myMaxSize != maxSize )
  {
    myMaxSize = maxSize;
    NotifySubMeshesHypothesisModification();
  }
}
//=======================================================================
//function : SaveTo
//purpose  : Persistence
ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save)
{
  save << myMinSize << " " << myMaxSize << " " << myDeflection;
  save << " " << -1 << " " << -1; // preview addition of parameters
  return save;
}
//=======================================================================
//function : LoadFrom
//purpose  : Persistence
istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
{
  int dummyParam;
  bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
  if (!isOK)
    load.clear(ios::badbit | load.rdstate());
  return load;
}
//=======================================================================
//function : SetParametersByMesh
//purpose  : Initialize parameters by the mesh built on the geometry
//param theMesh - the built mesh
//param theShape - the geometry of interest
//retval bool - true if parameter values have been successfully defined
bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh*   theMesh,
                                                const TopoDS_Shape& theShape)
{
  if ( !theMesh || theShape.IsNull() )
    return false;

  int nbEdges = 0;
  TopTools_IndexedMapOfShape edgeMap;
  TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );

  SMESH_MesherHelper helper( (SMESH_Mesh&) *theMesh );
  double minSz2 = 1e100, maxSz2 = 0, sz2, maxDefl2 = 0;
  for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
  {
    const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
    SMESHDS_SubMesh* smDS = theMesh->GetMeshDS()->MeshElements( edge );
    if ( !smDS ) continue;
    ++nbEdges;

    helper.SetSubShape( edge );
    BRepAdaptor_Curve curve( edge );

    SMDS_ElemIteratorPtr segIt = smDS->GetElements();
    while ( segIt->more() )
    {
      const SMDS_MeshElement* seg = segIt->next();
      const SMDS_MeshNode* n1 = seg->GetNode(0);
      const SMDS_MeshNode* n2 = seg->GetNode(1);
      sz2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
      minSz2 = Min( minSz2, sz2 );
      maxSz2 = Max( maxSz2, sz2 );
      if ( curve.GetType() != GeomAbs_Line )
      {
        double u1 = helper.GetNodeU( edge, n1, n2 );
        double u2 = helper.GetNodeU( edge, n2, n1 );
        maxDefl2 = Max( maxDefl2, deflection2( curve, u1, u2 ));
      }
    }
  }
  if ( nbEdges )
  {
    myMinSize = sqrt( minSz2 );
    myMaxSize = sqrt( maxSz2 );
    if ( maxDefl2 > 0 )
      myDeflection = maxDefl2;
  }
  return nbEdges;
}

//=======================================================================
//function : SetParametersByDefaults
//purpose  : Initialize my parameter values by default parameters.
//retval   : bool - true if parameter values have been successfully defined
bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults&  dflts,
                                                    const SMESH_Mesh* /*theMesh*/)
{
  myMinSize = dflts._elemLength / 100;
  myMaxSize = dflts._elemLength * 2;
  myDeflection = myMinSize / 10;
  return true;
}

//=======================================================================
//function : GetAlgo
//purpose  : Returns an algorithm that works using this hypothesis
//=======================================================================

StdMeshers_AdaptiveAlgo_1D* StdMeshers_Adaptive1D::GetAlgo() const
{
  if ( !myAlgo )
  {
    StdMeshers_AdaptiveAlgo_1D* newAlgo = 
      new StdMeshers_AdaptiveAlgo_1D( _gen->GetANewId(), _studyId, _gen );
    newAlgo->SetHypothesis( this );

    ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
  }
  return myAlgo;
}

//================================================================================
/*!
 * \brief Constructor
 */
//================================================================================

StdMeshers_AdaptiveAlgo_1D::StdMeshers_AdaptiveAlgo_1D(int        hypId,
                                                       int        studyId,
                                                       SMESH_Gen* gen)
  : StdMeshers_Regular_1D( hypId, studyId, gen ),
    myHyp(NULL)
{
  _name = "AdaptiveAlgo_1D";
}

//================================================================================
/*!
 * \brief Sets the hypothesis
 */
//================================================================================

void StdMeshers_AdaptiveAlgo_1D::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
{
  myHyp = hyp;
}

//================================================================================
/*!
 * \brief Creates segments on all given EDGEs
 */
//================================================================================

bool StdMeshers_AdaptiveAlgo_1D::Compute(SMESH_Mesh &         theMesh,
                                         const TopoDS_Shape & theShape,
                                         double*              theProgress,
                                         int*                 theProgressTic)
{
  *theProgress = 0.01;

  if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
    return error( "Bad parameters: min size > max size" );

  SMESH_MesherHelper helper( theMesh );
  const double grading = 0.7;

  TopTools_IndexedMapOfShape edgeMap, faceMap;
  TopExp::MapShapes( theShape,                 TopAbs_EDGE, edgeMap );
  TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );

  // Triangulate the shape with the given deflection ?????????
  Bnd_B3d box;
  {
    IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
    box = im.GetBox();
  }
  *theProgress = 0.3;

  // holder of segment size at each point
  SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );

  // minimal segment size that sizeTree can store with reasonable tree height
  const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );


  // working data of EDGEs
  vector< EdgeData > edges;
  {
    // sort EDGEs by length
    multimap< double, TopoDS_Edge > edgeOfLength;
    for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
    {
      const TopoDS_Edge & edge = TopoDS::Edge( edgeMap( iE ));
      if ( !SMESH_Algo::isDegenerated( edge) )
        edgeOfLength.insert( make_pair( EdgeLength( edge ), edge ));
    }
    edges.resize( edgeOfLength.size() );
    multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin();
    for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
    {
      const TopoDS_Edge & edge = len2edge->second;
      EdgeData& eData = edges[ iE ];
      eData.myC3d.Initialize( edge );
      eData.myLength  = EdgeLength( edge );
      eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
      eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
    }
  }

  // Take into account size of already existing segments
  SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
  while ( segIterator->more() )
  {
    const SMDS_MeshElement* seg = segIterator->next();
    sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
  }

  // Set size of segments according to the deflection

  StdMeshers_Regular_1D::_hypType = DEFLECTION;
  StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();

  list< double > params;
  for ( int iE = 0; iE < edges.size(); ++iE )
  {
    EdgeData& eData = edges[ iE ];
    //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;

    double f = eData.First().myU, l = eData.Last().myU;
    if ( !computeInternalParameters( theMesh, eData.myC3d, eData.myLength, f,l, params, false, false ))
      continue;
    if ( params.size() <= 1 && helper.IsClosedEdge( eData.Edge() ) ) // 2 segments on a circle
    {
      params.clear();
      for ( int i = 1; i < 6; ++i )
        params.push_back(( l - f ) * i/6. );
    }
    EdgeData::TPntIter where = --eData.myPoints.end();
    list< double >::const_iterator param = params.begin();
    for ( ; param != params.end(); ++param )
      eData.AddPoint( where, *param );

    EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
    for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
      sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
  }

  // Limit size of segments according to distance to closest FACE

  for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
  {
    const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
    // cout << "FACE " << iF << "/" << faceMap.Extent()
    //      << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;

    ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
    TriaTreeData*     triaSearcher = triaTree.GetTriaData();

    if ( BRepAdaptor_Surface( face ).GetType() != GeomAbs_Plane )
      triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );

    for ( int iE = 0; iE < edges.size(); ++iE )
    {
      EdgeData& eData = edges[ iE ];

      // check if the face is in topological contact with the edge
      bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
                         helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));

      bool sizeDecreased = true;
      for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable
      {
        // get points to check distance to the face
        EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
        pIt1->mySegSize = sizeTree.GetSize( pIt1->myP );
        for ( ; pIt2 != eData.myPoints.end(); )
        {
          pIt2->mySegSize = sizeTree.GetSize( pIt2->myP );
          double curSize = Min( pIt1->mySegSize, pIt2->mySegSize );
          if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
          {
            double midU  = 0.5*( pIt1->myU + pIt2->myU );
            gp_Pnt midP  = eData.myC3d.Value( midU );
            double midSz = sizeTree.GetSize( midP );
            pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
          }
          else
          {
            ++pIt1, ++pIt2;
          }
        }
        // check if the face is more distant than a half of the current segment size,
        // if not, segment size is decreased
        sizeDecreased = false;
        const gp_Pnt* avoidPnt = & eData.First().myP;
        for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end();  )
        {
          double distToFace =
            triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
          double allowedSize = Max( minSize, distToFace*( 1. + grading ));
          if ( 1.1 * allowedSize < pIt1->mySegSize  )
          {
            sizeDecreased = true;
            sizeTree.SetSize( pIt1->myP, allowedSize );
            // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
            //      << "\t SetSize " << allowedSize << " at "
            //      << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
            pIt2 = pIt1;
            if ( --pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
              sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
            pIt2 = pIt1;
            if ( ++pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
              sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
            pIt1->mySegSize = allowedSize;
          }
          ++pIt1;
          if ( & (*pIt1) == & eData.Last() )
            avoidPnt = & eData.Last().myP;
          else
            avoidPnt = NULL;

          if ( iLoop > 20 )
          {
#ifdef _DEBUG_
            cout << "Infinite loop in StdMeshers_AdaptiveAlgo_1D::Compute()" << endl;
#endif
            sizeDecreased = false;
            break;
          }
        }
      } // while ( sizeDecreased )
    } // loop on edges

    *theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
    if ( _computeCanceled )
      return false;

  } // loop on faceMap


  // Create segments

  SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
  _quadraticMesh = theMesh.GetHypothesis( edges[0].Edge(), quadHyp, /*andAncestors=*/true );
  helper.SetIsQuadratic( _quadraticMesh );

  for ( int iE = 0; iE < edges.size(); ++iE )
  {
    EdgeData& eData = edges[ iE ];

    // estimate roughly min segement size on the EDGE
    double edgeMinSize = myHyp->GetMaxSize();
    EdgeData::TPntIter pIt1 = eData.myPoints.begin();
    for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
      edgeMinSize = Min( edgeMinSize, sizeTree.GetSize( pIt1->myP ));

    const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
    const double parLen = l - f;
    const int  nbDivSeg = 5;
    int           nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );

    // compute nb of segments
    vector< double > nbSegs;
    bool toRecompute = true;
    double maxSegSize = 0;
    //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
    while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
    {
      nbSegs.resize( nbDiv + 1 );
      nbSegs[0] = 0;
      toRecompute = false;

      gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
      for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
      {
        p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
        double locSize = Min( sizeTree.GetSize( p2 ), myHyp->GetMaxSize() );
        double nb      = p1.Distance( p2 ) / locSize;
        // if ( nbSegs.size() < 30 )
        //   cout << "locSize " << locSize << " nb " << nb << endl;
        if ( nb > 1. )
        {
          toRecompute = true;
          edgeMinSize = locSize;
          nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
          break;
        }
        nbSegs[i] = nbSegs[i-1] + nb;
        p1 = p2;
        if ( nbSegs[i] >= segCount )
        {
          maxSegSize = Max( maxSegSize, pDiv.Distance( p2 ));
          pDiv = p2;
          ++segCount;
        }
      }
    }

    // compute parameters of nodes
    int nbSegFinal = int(floor(nbSegs.back()+0.5));
    double fact = nbSegFinal / nbSegs.back();
    if ( maxSegSize / fact > myHyp->GetMaxSize() )
      fact = ++nbSegFinal / nbSegs.back();
    //cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
    params.clear();
    for ( int i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
    {
      while ( nbSegs[i] * fact < segCount )
        ++i;
      if ( i < nbDiv )
        params.push_back( f + parLen * i / nbDiv );
      else
        break;
    }
    // get nodes on VERTEXes
    TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
    TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
    theMesh.GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
    theMesh.GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
    const SMDS_MeshNode * nf = VertexNode( vf, theMesh.GetMeshDS() );
    const SMDS_MeshNode * nl = VertexNode( vl, theMesh.GetMeshDS() );
    if ( !nf || !nl )
      return error("No node on vertex");

    // create segments
    helper.SetSubShape( eData.Edge() );
    helper.SetElementsOnShape( true );
    const int ID = 0;
    const SMDS_MeshNode *n1 = nf, *n2;
    list< double >::const_iterator u = params.begin();
    for ( ; u != params.end(); ++u, n1 = n2 )
    {
      gp_Pnt p2 = eData.myC3d.Value( *u );
      n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, *u );
      helper.AddEdge( n1, n2, ID, /*force3d=*/false );
    }
    helper.AddEdge( n1, nl, ID, /*force3d=*/false );

    *theProgress = 0.6 + 0.4 * iE / double( edges.size() );
    if ( _computeCanceled )
      return false;

  } // loop on EDGEs

}


//================================================================================
/*!
 * \brief Creates segments on all given EDGEs
 */
//================================================================================

bool StdMeshers_AdaptiveAlgo_1D::Evaluate(SMESH_Mesh &         theMesh,
                                          const TopoDS_Shape & theShape,
                                          MapShapeNbElems&     theResMap)
{
  // initialize fields of inherited StdMeshers_Regular_1D
  StdMeshers_Regular_1D::_hypType = DEFLECTION;
  StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();

  TopExp_Explorer edExp( theShape, TopAbs_EDGE );

  for ( ; edExp.More(); edExp.Next() )
  {
    const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
    StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );
  }
  return true;
}