src/StdMeshers/StdMeshers_Adaptive1D.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //  File   : StdMeshers_Adaptive1D.cxx
  23 //  Module : SMESH
  24 //
  25 #include "StdMeshers_Adaptive1D.hxx"
  26
  27 #include "SMESH_Gen.hxx"
  28 #include "SMESH_Mesh.hxx"
  29 #include "SMESH_MesherHelper.hxx"
  30 #include "SMESH_Octree.hxx"
  31 #include "SMESH_subMesh.hxx"
  32 #include "SMESH_HypoFilter.hxx"
  33
  34 #include <Utils_SALOME_Exception.hxx>
  35
  36 #include <BRepAdaptor_Curve.hxx>
  37 #include <BRepAdaptor_Surface.hxx>
  38 #include <BRepMesh_IncrementalMesh.hxx>
  39 #include <BRep_Tool.hxx>
  40 #include <Bnd_B3d.hxx>
  41 #include <GCPnts_AbscissaPoint.hxx>
  42 #include <GeomAdaptor_Curve.hxx>
  43 #include <Geom_Curve.hxx>
  44 #include <Poly_Array1OfTriangle.hxx>
  45 #include <Poly_Triangulation.hxx>
  46 #include <TColgp_Array1OfPnt.hxx>
  47 #include <TopExp.hxx>
  48 #include <TopExp_Explorer.hxx>
  49 #include <TopLoc_Location.hxx>
  50 #include <TopTools_IndexedMapOfShape.hxx>
  51 #include <TopoDS.hxx>
  52 #include <TopoDS_Edge.hxx>
  53 #include <TopoDS_Face.hxx>
  54 #include <TopoDS_Vertex.hxx>
  55 #include <gp_Lin.hxx>
  56 #include <gp_Pnt.hxx>
  57
  58 #include <limits>
  59 #include <vector>
  60 #include <set>
  61
  62 using namespace std;
  63
  64 namespace // internal utils
  65 {
  66   //================================================================================
  67   /*!
  68    * \brief Bnd_B3d with access to its center and half-size
  69    */
  70   struct BBox : public Bnd_B3d
  71   {
  72     gp_XYZ Center() const { return gp_XYZ( myCenter[0], myCenter[1], myCenter[2] ); }
  73     gp_XYZ HSize()  const { return gp_XYZ( myHSize[0],  myHSize[1],  myHSize[2]  ); }
  74     double Size()   const { return 2 * myHSize[0]; }
  75   };
  76   //================================================================================
  77   /*!
  78    * \brief Working data of an EDGE
  79    */
  80   struct EdgeData
  81   {
  82     struct ProbePnt
  83     {
  84       gp_Pnt myP;
  85       double myU;
  86       double mySegSize;
  87       ProbePnt( gp_Pnt p, double u, double sz=1e100 ): myP( p ), myU( u ), mySegSize( sz ) {}
  88     };
  89     BRepAdaptor_Curve myC3d;
  90     double            myLength;
  91     list< ProbePnt >  myPoints;
  92     BBox              myBBox;
  93
  94     typedef list< ProbePnt >::iterator TPntIter;
  95     void AddPoint( TPntIter where, double u )
  96     {
  97       TPntIter it = myPoints.insert( where, ProbePnt( myC3d.Value( u ), u ));
  98       myBBox.Add( it->myP.XYZ() );
  99     }
 100     const ProbePnt& First() const { return myPoints.front(); }
 101     const ProbePnt& Last()  const { return myPoints.back(); }
 102     const TopoDS_Edge& Edge() const { return myC3d.Edge(); }
 103     bool IsTooDistant( const BBox& faceBox, double maxSegSize ) const
 104     {
 105       gp_XYZ hsize = myBBox.HSize() + gp_XYZ( maxSegSize, maxSegSize, maxSegSize );
 106       return faceBox.IsOut ( Bnd_B3d( myBBox.Center(), hsize ));
 107     }
 108   };
 109   //================================================================================
 110   /*!
 111    * \brief Octree of local segment size
 112    */
 113   class SegSizeTree : public SMESH_Octree
 114   {
 115     double mySegSize; // segment size
 116
 117     // structure holding some common parameters of SegSizeTree
 118     struct _CommonData : public SMESH_TreeLimit
 119     {
 120       double myGrading, myMinSize, myMaxSize;
 121     };
 122     _CommonData* getData() const { return (_CommonData*) myLimit; }
 123
 124     SegSizeTree(double size): SMESH_Octree(), mySegSize(size)
 125     {
 126       allocateChildren();
 127     }
 128     void allocateChildren()
 129     {
 130       myChildren = new SMESH_Octree::TBaseTree*[nbChildren()];
 131       for ( int i = 0; i < nbChildren(); ++i )
 132         myChildren[i] = NULL;
 133     }
 134     virtual box_type* buildRootBox() { return 0; }
 135     virtual SegSizeTree* newChild() const { return 0; }
 136     virtual void buildChildrenData() {}
 137
 138   public:
 139
 140     SegSizeTree( Bnd_B3d & bb, double grading, double mixSize, double maxSize);
 141     void   SetSize( const gp_Pnt& p, double size );
 142     double SetSize( const gp_Pnt& p1, const gp_Pnt& p2 );
 143     double GetSize( const gp_Pnt& p ) const;
 144     const BBox* GetBox() const { return (BBox*) getBox(); }
 145     double GetMinSize() { return getData()->myMinSize; }
 146   };
 147   //================================================================================
 148   /*!
 149    * \brief Adaptive wire discertizator.
 150    */
 151   class AdaptiveAlgo : public StdMeshers_Regular_1D
 152   {
 153   public:
 154     AdaptiveAlgo(int hypId, int studyId, SMESH_Gen* gen);
 155     virtual bool Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape );
 156     virtual bool Evaluate(SMESH_Mesh &         theMesh,
 157                           const TopoDS_Shape & theShape,
 158                           MapShapeNbElems&     theResMap);
 159     void SetHypothesis( const StdMeshers_Adaptive1D* hyp );
 160   private:
 161
 162     bool makeSegments();
 163
 164     const StdMeshers_Adaptive1D* myHyp;
 165     SMESH_Mesh*                  myMesh;
 166     vector< EdgeData >           myEdges;
 167     SegSizeTree*                 mySizeTree;
 168   };
 169
 170   //================================================================================
 171   /*!
 172    * \brief Data of triangle used to locate it in an octree and to find distance
 173    *        to a point
 174    */
 175   struct Triangle
 176   {
 177     Bnd_B3d myBox;
 178     bool    myIsChecked; // to mark treated trias instead of using std::set
 179     // data for DistToProjection()
 180     gp_XYZ myN0, myEdge1, myEdge2, myNorm, myPVec;
 181     double myInvDet, myMaxSize2;
 182
 183     void Init( const gp_Pnt& n1, const gp_Pnt& n2, const gp_Pnt& n3 );
 184     bool DistToProjection( const gp_Pnt& p, double& dist ) const;
 185   };
 186   //================================================================================
 187   /*!
 188    * \brief Element data held by ElementBndBoxTree + algorithm computing a distance
 189    *        from a point to element
 190    */
 191   class ElementBndBoxTree;
 192   struct ElemTreeData : public SMESH_TreeLimit
 193   {
 194     vector< int >                myWorkIDs[8];// to speed up filling ElementBndBoxTree::_elementIDs
 195     virtual const Bnd_B3d* GetBox(int elemID) const = 0;
 196   };
 197   struct TriaTreeData : public ElemTreeData
 198   {
 199     vector< Triangle >           myTrias;
 200     double                       myFaceTol;
 201     double                       myTriasDeflection;
 202     BBox                         myBBox;
 203     BRepAdaptor_Surface          mySurface;
 204     ElementBndBoxTree*           myTree;
 205     const Poly_Array1OfTriangle* myPolyTrias;
 206     const TColgp_Array1OfPnt*    myNodes;
 207     bool                         myOwnNodes;
 208
 209     typedef vector<int> IntVec;
 210     IntVec                       myFoundTriaIDs;
 211
 212     TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree );
 213     ~TriaTreeData() { if ( myOwnNodes ) delete myNodes; myNodes = NULL; }
 214     virtual const Bnd_B3d* GetBox(int elemID) const { return &myTrias[elemID].myBox; }
 215     void PrepareToTriaSearch();
 216     void SetSizeByTrias( SegSizeTree& sizeTree, double deflection ) const;
 217     double GetMinDistInSphere(const gp_Pnt& p,
 218                               const double  radius,
 219                               const bool    projectedOnly,
 220                               const gp_Pnt* avoidP=0) const;
 221   };
 222   //================================================================================
 223   /*!
 224    * \brief Octree of triangles or segments
 225    */
 226   class ElementBndBoxTree : public SMESH_Octree
 227   {
 228   public:
 229     ElementBndBoxTree(const TopoDS_Face& face);
 230     void GetElementsInSphere( const gp_XYZ& center,
 231                               const double  radius, vector<int> & foundElemIDs) const;
 232     void FillIn();
 233     ElemTreeData* GetElemData() const { return (ElemTreeData*) myLimit; }
 234     TriaTreeData* GetTriaData() const { return (TriaTreeData*) myLimit; }
 235
 236   protected:
 237     ElementBndBoxTree() {}
 238     SMESH_Octree* newChild() const { return new ElementBndBoxTree; }
 239     void          buildChildrenData();
 240     Bnd_B3d*      buildRootBox();
 241   private:
 242     vector< int > _elementIDs;
 243   };
 244   //================================================================================
 245   /*!
 246    * \brief BRepMesh_IncrementalMesh with access to its protected Bnd_Box
 247    */
 248   struct IncrementalMesh : public BRepMesh_IncrementalMesh
 249   {
 250     IncrementalMesh(const TopoDS_Shape& shape,
 251                     const Standard_Real deflection,
 252                     const bool          relative):
 253       BRepMesh_IncrementalMesh( shape, deflection, relative )
 254     {
 255     }
 256     Bnd_B3d GetBox() const
 257     {
 258       Standard_Real TXmin, TYmin, TZmin, TXmax, TYmax, TZmax;
 259       myBox.Get(TXmin, TYmin, TZmin, TXmax, TYmax, TZmax);
 260       Bnd_B3d bb;
 261       bb.Add( gp_XYZ( TXmin, TYmin, TZmin ));
 262       bb.Add( gp_XYZ( TXmax, TYmax, TZmax ));
 263       return bb;
 264     }
 265   };
 266
 267   //================================================================================
 268   /*!
 269    * \brief Initialize TriaTreeData
 270    */
 271   //================================================================================
 272
 273   TriaTreeData::TriaTreeData( const TopoDS_Face& face, ElementBndBoxTree* triaTree )
 274     : myTriasDeflection(0), mySurface( face ),
 275       myTree(NULL), myPolyTrias(NULL), myNodes(NULL), myOwnNodes(false)
 276   {
 277     TopLoc_Location loc;
 278     Handle(Poly_Triangulation) tr = BRep_Tool::Triangulation( face, loc );
 279     if ( !tr.IsNull() )
 280     {
 281       myFaceTol         = SMESH_MesherHelper::MaxTolerance( face );
 282       myTree            = triaTree;
 283       myNodes           = & tr->Nodes();
 284       myPolyTrias       = & tr->Triangles();
 285       myTriasDeflection = tr->Deflection();
 286       if ( !loc.IsIdentity() ) // transform nodes if necessary
 287       {
 288         TColgp_Array1OfPnt* trsfNodes = new TColgp_Array1OfPnt( myNodes->Lower(), myNodes->Upper() );
 289         trsfNodes->Assign( *myNodes );
 290         myNodes    = trsfNodes;
 291         myOwnNodes = true;
 292         const gp_Trsf& trsf = loc;
 293         for ( int i = trsfNodes->Lower(); i <= trsfNodes->Upper(); ++i )
 294           trsfNodes->ChangeValue(i).Transform( trsf );
 295       }
 296       for ( int i = myNodes->Lower(); i <= myNodes->Upper(); ++i )
 297         myBBox.Add( myNodes->Value(i).XYZ() );
 298
 299     }
 300   }
 301   void TriaTreeData::PrepareToTriaSearch()
 302   {
 303     if ( !myTrias.empty() ) return; // already done
 304     if ( !myPolyTrias ) return;
 305
 306     myTrias.resize( myPolyTrias->Length() );
 307     Standard_Integer n1,n2,n3;
 308     for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
 309     {
 310       myPolyTrias->Value( i ).Get( n1,n2,n3 );
 311       myTrias[ i-1 ].Init( myNodes->Value( n1 ),
 312                            myNodes->Value( n2 ),
 313                            myNodes->Value( n3 ));
 314     }
 315     myTree->FillIn();
 316
 317     // TODO: mark triangles with nodes on VERTEXes to
 318     // less frequently compare with avoidPnt in GetMinDistInSphere()
 319     //
 320     // Handle(Poly_PolygonOnTriangulation) polygon =
 321     //   BRep_Tool::PolygonOnTriangulation( edge, tr, loc );
 322     // if ( polygon.IsNull() || !pologon.HasParameters() )
 323     //   continue;
 324     // Handle(TColStd_Array1OfInteger) nodeIDs = polygon->Nodes();
 325   }
 326
 327   //================================================================================
 328   /*!
 329    * \brief Set size of segments by size of triangles
 330    */
 331   //================================================================================
 332
 333   void TriaTreeData::SetSizeByTrias( SegSizeTree& sizeTree, double hypDeflection ) const
 334   {
 335     if ( mySurface.GetType() == GeomAbs_Plane ||
 336          myTriasDeflection   <= 1e-100 )
 337       return;
 338     const double factor = hypDeflection / myTriasDeflection;
 339
 340     bool isConstSize;
 341     switch( mySurface.GetType() ) {
 342     case GeomAbs_Cylinder:
 343     case GeomAbs_Sphere:
 344     case GeomAbs_Torus:
 345       isConstSize = true; break;
 346     default:
 347       isConstSize = false;
 348     }
 349
 350     typedef std::pair<int,int> TLink;
 351     TLink link;
 352     map< TLink, double >           lenOfDoneLink;
 353     map< TLink, double >::iterator link2len;
 354
 355     Standard_Integer n[4];
 356     gp_Pnt p[4];
 357     double a[3];
 358     bool   isDone[3];
 359     double size = -1., maxLinkLen;
 360     int    jLongest;
 361
 362     //int nbLinks = 0;
 363     for ( int i = 1; i <= myPolyTrias->Upper(); ++i )
 364     {
 365       // get corners of a triangle
 366       myPolyTrias->Value( i ).Get( n[0],n[1],n[2] );
 367       n[3] = n[0];
 368       p[0] = myNodes->Value( n[0] );
 369       p[1] = myNodes->Value( n[1] );
 370       p[2] = myNodes->Value( n[2] );
 371       p[3] = p[0];
 372       // get length of links and find the longest one
 373       maxLinkLen = 0;
 374       for ( int j = 0; j < 3; ++j )
 375       {
 376         if ( n[j] < n[j+1] )
 377           link = TLink( n[j], n[j+1] );
 378         else
 379           link = TLink( n[j+1], n[j] );
 380         link2len  = lenOfDoneLink.insert( make_pair( link, -1. )).first;
 381         isDone[j] = !((*link2len).second < 0 );
 382         a[j]      = isDone[j] ? (*link2len).second : (*link2len).second = p[j].Distance( p[j+1] );
 383         if ( isDone[j] )
 384           lenOfDoneLink.erase( link2len );
 385         if ( a[j] > maxLinkLen )
 386         {
 387           maxLinkLen = a[j];
 388           jLongest   = j;
 389         }
 390       }
 391       // compute minimal altitude of a triangle
 392       if ( !isConstSize || size < 0. )
 393       {
 394         double s    = 0.5 * ( a[0] + a[1] + a[2] );
 395         double area = sqrt( s * (s - a[0]) * (s - a[1]) * (s - a[2]));
 396         size        = 2 * area / maxLinkLen; // minimal altitude
 397       }
 398       // set size to the size tree
 399       if ( !isDone[ jLongest ] || !isConstSize )
 400       {
 401         //++nbLinks;
 402         if ( size < numeric_limits<double>::min() )
 403           continue;
 404         int nb = Max( 1, int( maxLinkLen / size / 2 ));
 405         for ( int k = 0; k <= nb; ++k )
 406         {
 407           double r = double( k ) / nb;
 408           sizeTree.SetSize( r * p[ jLongest ].XYZ() + ( 1-r ) * p[ jLongest+1 ].XYZ(),
 409                             size * factor );
 410         }
 411       }
 412       //cout << "SetSizeByTrias, i="<< i << " " << sz * factor << endl;
 413     }
 414     // cout << "SetSizeByTrias, nn tria="<< myPolyTrias->Upper()
 415     //      << " nb links" << nbLinks << " isConstSize="<<isConstSize
 416     //      << " " << size * factor << endl;
 417   }
 418   //================================================================================
 419   /*!
 420    * \brief Return minimal distance from a given point to a trinangle but not more
 421    *        distant than a given radius. Triangles with a node at avoidPnt are ignored.
 422    *        If projectedOnly,
 423    */
 424   //================================================================================
 425
 426   double TriaTreeData::GetMinDistInSphere(const gp_Pnt& p,
 427                                           const double  radius,
 428                                           const bool    projectedOnly,
 429                                           const gp_Pnt* avoidPnt) const
 430   {
 431     double minDist2 = 1e100;
 432     const double tol2 = myFaceTol * myFaceTol;
 433
 434     TriaTreeData* me = const_cast<TriaTreeData*>( this );
 435     me->myFoundTriaIDs.clear();
 436     myTree->GetElementsInSphere( p.XYZ(), radius, me->myFoundTriaIDs );
 437
 438     Standard_Integer n[ 3 ];
 439     for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
 440     {
 441       Triangle& t = me->myTrias[ myFoundTriaIDs[i] ];
 442       if ( t.myIsChecked )
 443         continue;
 444       t.myIsChecked = true;
 445
 446       double d, minD2 = minDist2;
 447       bool avoidTria = false;
 448       myPolyTrias->Value( myFoundTriaIDs[i]+1 ).Get( n[0],n[1],n[2] );
 449       for ( int i = 0; i < 3; ++i )
 450       {
 451         const gp_Pnt& pn = myNodes->Value(n[i]);
 452         if ( avoidTria = ( avoidPnt && pn.SquareDistance(*avoidPnt) <= tol2 ))
 453           break;
 454         if ( !projectedOnly )
 455           minD2 = Min( minD2, pn.SquareDistance( p ));
 456       }
 457       if ( !avoidTria )
 458       {
 459         if ( minD2 < t.myMaxSize2 && t.DistToProjection( p, d ))
 460           minD2 = Min( minD2, d*d );
 461         minDist2 = Min( minDist2, minD2 );
 462       }
 463     }
 464
 465     for ( size_t i = 0; i < myFoundTriaIDs.size(); ++i )
 466       me->myTrias[ myFoundTriaIDs[i] ].myIsChecked = false;
 467
 468     return sqrt( minDist2 );
 469   }
 470   //================================================================================
 471   /*!
 472    * \brief Prepare Triangle data
 473    */
 474   //================================================================================
 475
 476   void Triangle::Init( const gp_Pnt& p1, const gp_Pnt& p2, const gp_Pnt& p3 )
 477   {
 478     myBox.Add( p1 );
 479     myBox.Add( p2 );
 480     myBox.Add( p3 );
 481     myN0 = p1.XYZ();
 482     myEdge1 = p2.XYZ() - myN0;
 483     myEdge2 = p3.XYZ() - myN0;
 484     myNorm = myEdge1 ^ myEdge2;
 485     double normSize = myNorm.Modulus();
 486     if ( normSize > std::numeric_limits<double>::min() )
 487     {
 488       myNorm /= normSize;
 489       myPVec = myNorm ^ myEdge2;
 490       myInvDet = 1. / ( myEdge1 * myPVec );
 491     }
 492     else
 493     {
 494       myInvDet = 0.;
 495     }
 496     myMaxSize2 = Max( p2.SquareDistance( p3 ),
 497                       Max( myEdge2.SquareModulus(), myEdge1.SquareModulus() ));
 498   }
 499   //================================================================================
 500   /*!
 501    * \brief Compute distance from a point to the triangle. Return false if the point
 502    *        is not projected inside the triangle
 503    */
 504   //================================================================================
 505
 506   bool Triangle::DistToProjection( const gp_Pnt& p, double& dist ) const
 507   {
 508     if ( myInvDet == 0 )
 509       return false; // degenerated triangle
 510
 511     /* distance from n0 to the point */
 512     gp_XYZ tvec = p.XYZ() - myN0;
 513
 514     /* calculate U parameter and test bounds */
 515     double u = ( tvec * myPVec ) * myInvDet;
 516     if (u < 0.0 || u > 1.0)
 517       return false; // projected outside the triangle
 518
 519     /* calculate V parameter and test bounds */
 520     gp_XYZ qvec = tvec ^ myEdge1;
 521     double v = ( myNorm * qvec) * myInvDet;
 522     if ( v < 0.0 || u + v > 1.0 )
 523       return false; // projected outside the triangle
 524
 525     dist = ( myEdge2 * qvec ) * myInvDet;
 526     return true;
 527   }
 528
 529   //================================================================================
 530   /*!
 531    * \brief Consturct ElementBndBoxTree of Poly_Triangulation of a FACE
 532    */
 533   //================================================================================
 534
 535   ElementBndBoxTree::ElementBndBoxTree(const TopoDS_Face& face)
 536     :SMESH_Octree()
 537   {
 538     TriaTreeData* data = new TriaTreeData( face, this );
 539     data->myMaxLevel = 5;
 540     myLimit = data;
 541   }
 542   //================================================================================
 543   /*!
 544    * \brief Fill all levels of octree of Poly_Triangulation of a FACE
 545    */
 546   //================================================================================
 547
 548   void ElementBndBoxTree::FillIn()
 549   {
 550     if ( myChildren ) return;
 551     TriaTreeData* data = GetTriaData();
 552     if ( !data->myTrias.empty() )
 553     {
 554       for ( size_t i = 0; i < data->myTrias.size(); ++i )
 555         _elementIDs.push_back( i );
 556
 557       compute();
 558     }
 559   }
 560   //================================================================================
 561   /*!
 562    * \brief Return the maximal box
 563    */
 564   //================================================================================
 565
 566   Bnd_B3d* ElementBndBoxTree::buildRootBox()
 567   {
 568     TriaTreeData* data = GetTriaData();
 569     Bnd_B3d*       box = new Bnd_B3d( data->myBBox );
 570     return box;
 571   }
 572   //================================================================================
 573   /*!
 574    * \brief Redistrubute element boxes among children
 575    */
 576   //================================================================================
 577
 578   void ElementBndBoxTree::buildChildrenData()
 579   {
 580     ElemTreeData* data = GetElemData();
 581     for ( int i = 0; i < _elementIDs.size(); ++i )
 582     {
 583       const Bnd_B3d* elemBox = data->GetBox( _elementIDs[i] );
 584       for (int j = 0; j < 8; j++)
 585         if ( !elemBox->IsOut( *myChildren[ j ]->getBox() ))
 586           data->myWorkIDs[ j ].push_back( _elementIDs[i] );
 587     }
 588     SMESHUtils::FreeVector( _elementIDs ); // = _elements.clear() + free memory
 589
 590     const int theMaxNbElemsInLeaf = 7;
 591
 592     for (int j = 0; j < 8; j++)
 593     {
 594       ElementBndBoxTree* child = static_cast<ElementBndBoxTree*>( myChildren[j] );
 595       child->_elementIDs = data->myWorkIDs[ j ];
 596       if ( child->_elementIDs.size() <= theMaxNbElemsInLeaf )
 597         child->myIsLeaf = true;
 598       data->myWorkIDs[ j ].clear();
 599     }
 600   }
 601   //================================================================================
 602   /*!
 603    * \brief Return elements from leaves intersecting the sphere
 604    */
 605   //================================================================================
 606
 607   void ElementBndBoxTree::GetElementsInSphere( const gp_XYZ& center,
 608                                                const double  radius,
 609                                                vector<int> & foundElemIDs) const
 610   {
 611     if ( const box_type* box = getBox() )
 612     {
 613       if ( box->IsOut( center, radius ))
 614         return;
 615
 616       if ( isLeaf() )
 617       {
 618         ElemTreeData* data = GetElemData();
 619         for ( int i = 0; i < _elementIDs.size(); ++i )
 620           if ( !data->GetBox( _elementIDs[i] )->IsOut( center, radius ))
 621             foundElemIDs.push_back( _elementIDs[i] );
 622       }
 623       else
 624       {
 625         for (int i = 0; i < 8; i++)
 626           ((ElementBndBoxTree*) myChildren[i])->GetElementsInSphere( center, radius, foundElemIDs );
 627       }
 628     }
 629   }
 630
 631   //================================================================================
 632   /*!
 633    * \brief Constructor of SegSizeTree
 634    *  \param [in,out] bb - bounding box enclosing all EDGEs to discretize
 635    *  \param [in] grading - factor to get max size of the neighbour segment by
 636    *         size of a current one.
 637    */
 638   //================================================================================
 639
 640   SegSizeTree::SegSizeTree( Bnd_B3d & bb, double grading, double minSize, double maxSize )
 641     : SMESH_Octree( new _CommonData() )
 642   {
 643     // make cube myBox from the box bb
 644     gp_XYZ pmin = bb.CornerMin(), pmax = bb.CornerMax();
 645     double maxBoxHSize = 0.5 * Max( pmax.X()-pmin.X(), Max( pmax.Y()-pmin.Y(), pmax.Z()-pmin.Z() ));
 646     maxBoxHSize *= 1.01;
 647     bb.SetHSize( gp_XYZ( maxBoxHSize, maxBoxHSize, maxBoxHSize ));
 648     myBox = new box_type( bb );
 649
 650     mySegSize = Min( 2 * maxBoxHSize, maxSize );
 651
 652     getData()->myGrading = grading;
 653     getData()->myMinSize = Max( minSize, 2*maxBoxHSize / 1.e6 );
 654     getData()->myMaxSize = maxSize;
 655     allocateChildren();
 656   }
 657
 658   //================================================================================
 659   /*!
 660    * \brief Set segment size at a given point
 661    */
 662   //================================================================================
 663
 664   void SegSizeTree::SetSize( const gp_Pnt& p, double size )
 665   {
 666     // check if the point is out of the largest cube
 667     SegSizeTree* root = this;
 668     while ( root->myFather )
 669       root = (SegSizeTree*) root->myFather;
 670     if ( root->getBox()->IsOut( p.XYZ() ))
 671       return;
 672
 673     // keep size whthin the valid range
 674     size = Max( size, getData()->myMinSize );
 675     //size = Min( size, getData()->myMaxSize );
 676
 677     // find an existing leaf at the point
 678     SegSizeTree* leaf = (SegSizeTree*) root;
 679     int iChild;
 680     while ( true )
 681     {
 682       iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
 683       if ( leaf->myChildren[ iChild ] )
 684         leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
 685       else
 686         break;
 687     }
 688     // don't increase the current size
 689     if ( leaf->mySegSize <= 1.1 * size )
 690       return;
 691
 692     // split the found leaf until its box size is less than the given size
 693     const double rootSize = root->GetBox()->Size();
 694     while ( leaf->GetBox()->Size() > size )
 695     {
 696       const BBox* bb = leaf->GetBox();
 697       iChild   = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), bb->Center() );
 698       SegSizeTree* newLeaf = new SegSizeTree( bb->Size() / 2 );
 699       leaf->myChildren[iChild] = newLeaf;
 700       newLeaf->myFather = leaf;
 701       newLeaf->myLimit  = leaf->myLimit;
 702       newLeaf->myLevel  = leaf->myLevel + 1;
 703       newLeaf->myBox    = leaf->newChildBox( iChild );
 704       newLeaf->myBox->Enlarge( rootSize * 1e-10 );
 705       //newLeaf->myIsLeaf = ( newLeaf->mySegSize <= size );
 706       leaf = newLeaf;
 707     }
 708     leaf->mySegSize = size;
 709
 710     // propagate increased size out from the leaf
 711     double boxSize = leaf->GetBox()->Size();
 712     double sizeInc = size + boxSize * getData()->myGrading;
 713     for ( int iDir = 1; iDir <= 3; ++iDir )
 714     {
 715       gp_Pnt outPnt = p;
 716       outPnt.SetCoord( iDir, p.Coord( iDir ) + boxSize );
 717       SetSize( outPnt, sizeInc );
 718       outPnt.SetCoord( iDir, p.Coord( iDir ) - boxSize );
 719       SetSize( outPnt, sizeInc );
 720     }
 721   }
 722   //================================================================================
 723   /*!
 724    * \brief Set size of a segment given by two end points
 725    */
 726   //================================================================================
 727
 728   double SegSizeTree::SetSize( const gp_Pnt& p1, const gp_Pnt& p2 )
 729   {
 730     const double size = p1.Distance( p2 );
 731     gp_XYZ p = 0.5 * ( p1.XYZ() + p2.XYZ() );
 732     SetSize( p, size );
 733     SetSize( p1, size );
 734     SetSize( p2, size );
 735     //cout << "SetSize " << p1.Distance( p2 ) << " at " << p.X() <<", "<< p.Y()<<", "<<p.Z()<< endl;
 736     return GetSize( p );
 737   }
 738
 739   //================================================================================
 740   /*!
 741    * \brief Return segment size at a point
 742    */
 743   //================================================================================
 744
 745   double SegSizeTree::GetSize( const gp_Pnt& p ) const
 746   {
 747     const SegSizeTree* leaf = this;
 748     while ( true )
 749     {
 750       int iChild = SMESH_Octree::getChildIndex( p.X(), p.Y(), p.Z(), leaf->GetBox()->Center() );
 751       if ( leaf->myChildren[ iChild ] )
 752         leaf = (SegSizeTree*) leaf->myChildren[ iChild ];
 753       else
 754         return leaf->mySegSize;
 755     }
 756     return mySegSize; // just to return anything
 757   }
 758
 759   //================================================================================
 760   /*!
 761    * \brief Evaluate curve deflection between two points
 762    * \param theCurve - the curve
 763    * \param theU1 - the parameter of the first point
 764    * \param theU2 - the parameter of the second point
 765    * \retval double - square deflection value
 766    */
 767   //================================================================================
 768
 769   double deflection2(const BRepAdaptor_Curve & theCurve,
 770                      double                    theU1,
 771                      double                    theU2)
 772   {
 773     // line between theU1 and theU2
 774     gp_Pnt p1 = theCurve.Value( theU1 ), p2 = theCurve.Value( theU2 );
 775     gp_Lin segment( p1, gp_Vec( p1, p2 ));
 776
 777     // evaluate square distance of theCurve from the segment
 778     Standard_Real dist2 = 0;
 779     const int nbPnt = 5;
 780     const double step = ( theU2 - theU1 ) / nbPnt;
 781     while (( theU1 += step ) < theU2 )
 782       dist2 = Max( dist2, segment.SquareDistance( theCurve.Value( theU1 )));
 783
 784     return dist2;
 785   }
 786
 787 } // namespace
 788
 789 //=======================================================================
 790 //function : StdMeshers_Adaptive1D
 791 //purpose  : Constructor
 792 StdMeshers_Adaptive1D::StdMeshers_Adaptive1D(int         hypId,
 793                                              int         studyId,
 794                                              SMESH_Gen * gen)
 795   :SMESH_Hypothesis(hypId, studyId, gen)
 796 {
 797   myMinSize       = 1e-10;
 798   myMaxSize       = 1e+10;
 799   myDeflection    = 1e-2;
 800   myAlgo          = NULL;
 801   _name           = "Adaptive1D";
 802   _param_algo_dim = 1; // is used by SMESH_Regular_1D
 803 }
 804 //=======================================================================
 805 //function : ~StdMeshers_Adaptive1D
 806 //purpose  : Destructor
 807 StdMeshers_Adaptive1D::~StdMeshers_Adaptive1D()
 808 {
 809   delete myAlgo; myAlgo = NULL;
 810 }
 811 //=======================================================================
 812 //function : SetDeflection
 813 //purpose  :
 814 void StdMeshers_Adaptive1D::SetDeflection(double value)
 815   throw(SALOME_Exception)
 816 {
 817   if (value <= std::numeric_limits<double>::min() )
 818     throw SALOME_Exception("Deflection must be greater that zero");
 819   if (myDeflection != value)
 820   {
 821     myDeflection = value;
 822     NotifySubMeshesHypothesisModification();
 823   }
 824 }
 825 //=======================================================================
 826 //function : SetMinSize
 827 //purpose  : Sets minimal allowed segment length
 828 void StdMeshers_Adaptive1D::SetMinSize(double minSize)
 829   throw(SALOME_Exception)
 830 {
 831   if (minSize <= std::numeric_limits<double>::min() )
 832     throw SALOME_Exception("Min size must be greater that zero");
 833
 834   if (myMinSize != minSize )
 835   {
 836     myMinSize = minSize;
 837     NotifySubMeshesHypothesisModification();
 838   }
 839 }
 840 //=======================================================================
 841 //function : SetMaxSize
 842 //purpose  : Sets maximal allowed segment length
 843 void StdMeshers_Adaptive1D::SetMaxSize(double maxSize)
 844   throw(SALOME_Exception)
 845 {
 846   if (maxSize <= std::numeric_limits<double>::min() )
 847     throw SALOME_Exception("Max size must be greater that zero");
 848
 849   if (myMaxSize != maxSize )
 850   {
 851     myMaxSize = maxSize;
 852     NotifySubMeshesHypothesisModification();
 853   }
 854 }
 855 //=======================================================================
 856 //function : SaveTo
 857 //purpose  : Persistence
 858 ostream & StdMeshers_Adaptive1D::SaveTo(ostream & save)
 859 {
 860   save << myMinSize << " " << myMaxSize << " " << myDeflection;
 861   save << " " << -1 << " " << -1; // preview addition of parameters
 862   return save;
 863 }
 864 //=======================================================================
 865 //function : LoadFrom
 866 //purpose  : Persistence
 867 istream & StdMeshers_Adaptive1D::LoadFrom(istream & load)
 868 {
 869   int dummyParam;
 870   bool isOK = (load >> myMinSize >> myMaxSize >> myDeflection >> dummyParam >> dummyParam);
 871   if (!isOK)
 872     load.clear(ios::badbit | load.rdstate());
 873   return load;
 874 }
 875 //=======================================================================
 876 //function : SetParametersByMesh
 877 //purpose  : Initialize parameters by the mesh built on the geometry
 878 //param theMesh - the built mesh
 879 //param theShape - the geometry of interest
 880 //retval bool - true if parameter values have been successfully defined
 881 bool StdMeshers_Adaptive1D::SetParametersByMesh(const SMESH_Mesh*   theMesh,
 882                                                 const TopoDS_Shape& theShape)
 883 {
 884   if ( !theMesh || theShape.IsNull() )
 885     return false;
 886
 887   int nbEdges = 0;
 888   TopTools_IndexedMapOfShape edgeMap;
 889   TopExp::MapShapes( theShape, TopAbs_EDGE, edgeMap );
 890
 891   SMESH_MesherHelper helper( (SMESH_Mesh&) *theMesh );
 892   double minSz2 = 1e100, maxSz2 = 0, sz2, maxDefl2 = 0;
 893   for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
 894   {
 895     const TopoDS_Edge& edge = TopoDS::Edge( edgeMap( iE ));
 896     SMESHDS_SubMesh* smDS = theMesh->GetMeshDS()->MeshElements( edge );
 897     if ( !smDS ) continue;
 898     ++nbEdges;
 899
 900     helper.SetSubShape( edge );
 901     BRepAdaptor_Curve curve( edge );
 902
 903     SMDS_ElemIteratorPtr segIt = smDS->GetElements();
 904     while ( segIt->more() )
 905     {
 906       const SMDS_MeshElement* seg = segIt->next();
 907       const SMDS_MeshNode* n1 = seg->GetNode(0);
 908       const SMDS_MeshNode* n2 = seg->GetNode(1);
 909       sz2 = SMESH_TNodeXYZ( n1 ).SquareDistance( n2 );
 910       minSz2 = Min( minSz2, sz2 );
 911       maxSz2 = Max( maxSz2, sz2 );
 912       if ( curve.GetType() != GeomAbs_Line )
 913       {
 914         double u1 = helper.GetNodeU( edge, n1, n2 );
 915         double u2 = helper.GetNodeU( edge, n2, n1 );
 916         maxDefl2 = Max( maxDefl2, deflection2( curve, u1, u2 ));
 917       }
 918     }
 919   }
 920   if ( nbEdges )
 921   {
 922     myMinSize = sqrt( minSz2 );
 923     myMaxSize = sqrt( maxSz2 );
 924     if ( maxDefl2 > 0 )
 925       myDeflection = maxDefl2;
 926   }
 927   return nbEdges;
 928 }
 929
 930 //=======================================================================
 931 //function : SetParametersByDefaults
 932 //purpose  : Initialize my parameter values by default parameters.
 933 //retval   : bool - true if parameter values have been successfully defined
 934 bool StdMeshers_Adaptive1D::SetParametersByDefaults(const TDefaults&  dflts,
 935                                                     const SMESH_Mesh* /*theMesh*/)
 936 {
 937   myMinSize = dflts._elemLength / 10;
 938   myMaxSize = dflts._elemLength * 2;
 939   myDeflection = myMinSize / 7;
 940   return true;
 941 }
 942
 943 //=======================================================================
 944 //function : GetAlgo
 945 //purpose  : Returns an algorithm that works using this hypothesis
 946 //=======================================================================
 947
 948 SMESH_Algo* StdMeshers_Adaptive1D::GetAlgo() const
 949 {
 950   if ( !myAlgo )
 951   {
 952     AdaptiveAlgo* newAlgo =
 953       new AdaptiveAlgo( _gen->GetANewId(), _studyId, _gen );
 954     newAlgo->SetHypothesis( this );
 955
 956     ((StdMeshers_Adaptive1D*) this)->myAlgo = newAlgo;
 957   }
 958   return myAlgo;
 959 }
 960
 961 //================================================================================
 962 /*!
 963  * \brief Constructor
 964  */
 965 //================================================================================
 966
 967 AdaptiveAlgo::AdaptiveAlgo(int        hypId,
 968                            int        studyId,
 969                            SMESH_Gen* gen)
 970   : StdMeshers_Regular_1D( hypId, studyId, gen ),
 971     myHyp(NULL)
 972 {
 973   _name = "AdaptiveAlgo_1D";
 974 }
 975
 976 //================================================================================
 977 /*!
 978  * \brief Sets the hypothesis
 979  */
 980 //================================================================================
 981
 982 void AdaptiveAlgo::SetHypothesis( const StdMeshers_Adaptive1D* hyp )
 983 {
 984   myHyp = hyp;
 985 }
 986
 987 //================================================================================
 988 /*!
 989  * \brief Creates segments on all given EDGEs
 990  */
 991 //================================================================================
 992
 993 bool AdaptiveAlgo::Compute(SMESH_Mesh &         theMesh,
 994                            const TopoDS_Shape & theShape)
 995 {
 996   //*theProgress = 0.01;
 997
 998   if ( myHyp->GetMinSize() > myHyp->GetMaxSize() )
 999     return error( "Bad parameters: min size > max size" );
1000
1001   myMesh = &theMesh;
1002   SMESH_MesherHelper helper( theMesh );
1003   const double grading = 0.7;
1004
1005   TopTools_IndexedMapOfShape edgeMap, faceMap;
1006   TopExp::MapShapes( theShape,                 TopAbs_EDGE, edgeMap );
1007   TopExp::MapShapes( theMesh.GetShapeToMesh(), TopAbs_FACE, faceMap );
1008
1009   // Triangulate the shape with the given deflection ?????????
1010   Bnd_B3d box;
1011   {
1012     IncrementalMesh im( theMesh.GetShapeToMesh(), myHyp->GetDeflection(), /*Relatif=*/false);
1013     box = im.GetBox();
1014   }
1015   //*theProgress = 0.3;
1016
1017   // holder of segment size at each point
1018   SegSizeTree sizeTree( box, grading, myHyp->GetMinSize(), myHyp->GetMaxSize() );
1019   mySizeTree = & sizeTree;
1020
1021   // minimal segment size that sizeTree can store with reasonable tree height
1022   const double minSize = Max( myHyp->GetMinSize(), 1.1 * sizeTree.GetMinSize() );
1023
1024
1025   // fill myEdges - working data of EDGEs
1026   {
1027     // sort EDGEs by length
1028     multimap< double, TopoDS_Edge > edgeOfLength;
1029     for ( int iE = 1; iE <= edgeMap.Extent(); ++iE )
1030     {
1031       const TopoDS_Edge & edge = TopoDS::Edge( edgeMap( iE ));
1032       if ( !SMESH_Algo::isDegenerated( edge) )
1033         edgeOfLength.insert( make_pair( EdgeLength( edge ), edge ));
1034     }
1035     myEdges.clear();
1036     myEdges.resize( edgeOfLength.size() );
1037     multimap< double, TopoDS_Edge >::const_iterator len2edge = edgeOfLength.begin();
1038     for ( int iE = 0; len2edge != edgeOfLength.end(); ++len2edge, ++iE )
1039     {
1040       const TopoDS_Edge & edge = len2edge->second;
1041       EdgeData& eData = myEdges[ iE ];
1042       eData.myC3d.Initialize( edge );
1043       eData.myLength  = EdgeLength( edge );
1044       eData.AddPoint( eData.myPoints.end(), eData.myC3d.FirstParameter() );
1045       eData.AddPoint( eData.myPoints.end(), eData.myC3d.LastParameter() );
1046     }
1047   }
1048   if ( _computeCanceled ) return false;
1049
1050   // Take into account size of already existing segments
1051   SMDS_EdgeIteratorPtr segIterator = theMesh.GetMeshDS()->edgesIterator();
1052   while ( segIterator->more() )
1053   {
1054     const SMDS_MeshElement* seg = segIterator->next();
1055     sizeTree.SetSize( SMESH_TNodeXYZ( seg->GetNode( 0 )), SMESH_TNodeXYZ( seg->GetNode( 1 )));
1056   }
1057   if ( _computeCanceled ) return false;
1058
1059   // Set size of segments according to the deflection
1060
1061   StdMeshers_Regular_1D::_hypType = DEFLECTION;
1062   StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1063
1064   list< double > params;
1065   for ( int iE = 0; iE < myEdges.size(); ++iE )
1066   {
1067     EdgeData& eData = myEdges[ iE ];
1068     //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1069
1070     double f = eData.First().myU, l = eData.Last().myU;
1071     if ( !computeInternalParameters( theMesh, eData.myC3d, eData.myLength, f,l, params, false, false ))
1072       continue;
1073     if ( params.size() <= 1 && helper.IsClosedEdge( eData.Edge() ) ) // 2 segments on a circle
1074     {
1075       params.clear();
1076       for ( int i = 1; i < 6; ++i )
1077         params.push_back(( l - f ) * i/6. );
1078     }
1079     EdgeData::TPntIter where = --eData.myPoints.end();
1080     list< double >::const_iterator param = params.begin();
1081     for ( ; param != params.end(); ++param )
1082       eData.AddPoint( where, *param );
1083
1084     EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1085     for ( ; pIt2 != eData.myPoints.end(); ++pIt1, ++pIt2 )
1086     {
1087       double sz = sizeTree.SetSize( (*pIt1).myP, (*pIt2).myP );
1088       sz = Min( sz, myHyp->GetMaxSize() );
1089       pIt1->mySegSize = Min( sz, pIt1->mySegSize );
1090       pIt2->mySegSize = Min( sz, pIt2->mySegSize );
1091     }
1092
1093     if ( _computeCanceled ) return false;
1094   }
1095
1096   // Limit size of segments according to distance to closest FACE
1097
1098   for ( int iF = 1; iF <= faceMap.Extent(); ++iF )
1099   {
1100     if ( _computeCanceled ) return false;
1101
1102     const TopoDS_Face & face = TopoDS::Face( faceMap( iF ));
1103     // cout << "FACE " << iF << "/" << faceMap.Extent()
1104     //      << " id-" << theMesh.GetMeshDS()->ShapeToIndex( face ) << endl;
1105
1106     ElementBndBoxTree triaTree( face ); // tree of FACE triangulation
1107     TriaTreeData*     triaSearcher = triaTree.GetTriaData();
1108
1109     triaSearcher->SetSizeByTrias( sizeTree, myHyp->GetDeflection() );
1110
1111     for ( int iE = 0; iE < myEdges.size(); ++iE )
1112     {
1113       EdgeData& eData = myEdges[ iE ];
1114
1115       // check if the face is in topological contact with the edge
1116       bool isAdjFace = ( helper.IsSubShape( helper.IthVertex( 0, eData.Edge()), face ) ||
1117                          helper.IsSubShape( helper.IthVertex( 1, eData.Edge()), face ));
1118
1119       if ( isAdjFace && triaSearcher->mySurface.GetType() == GeomAbs_Plane )
1120         continue;
1121
1122       bool sizeDecreased = true;
1123       for (int iLoop = 0; sizeDecreased; ++iLoop ) //repeat until segment size along the edge becomes stable
1124       {
1125         double maxSegSize = 0;
1126
1127         // get points to check distance to the face
1128         EdgeData::TPntIter pIt2 = eData.myPoints.begin(), pIt1 = pIt2++;
1129         maxSegSize = pIt1->mySegSize = Min( pIt1->mySegSize, sizeTree.GetSize( pIt1->myP ));
1130         for ( ; pIt2 != eData.myPoints.end(); )
1131         {
1132           pIt2->mySegSize = Min( pIt2->mySegSize, sizeTree.GetSize( pIt2->myP ));
1133           double curSize  = Min( pIt1->mySegSize, pIt2->mySegSize );
1134           maxSegSize      = Max( pIt2->mySegSize, maxSegSize );
1135           if ( pIt1->myP.Distance( pIt2->myP ) > curSize )
1136           {
1137             double midU  = 0.5*( pIt1->myU + pIt2->myU );
1138             gp_Pnt midP  = eData.myC3d.Value( midU );
1139             double midSz = sizeTree.GetSize( midP );
1140             pIt2 = eData.myPoints.insert( pIt2, EdgeData::ProbePnt( midP, midU, midSz ));
1141             eData.myBBox.Add( midP.XYZ() );
1142           }
1143           else
1144           {
1145             ++pIt1, ++pIt2;
1146           }
1147         }
1148         // check if the face is more distant than a half of the current segment size,
1149         // if not, segment size is decreased
1150
1151         if ( iLoop == 0 && eData.IsTooDistant( triaSearcher->myBBox, maxSegSize ))
1152           break;
1153         triaSearcher->PrepareToTriaSearch();
1154
1155         //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1156         sizeDecreased = false;
1157         const gp_Pnt* avoidPnt = & eData.First().myP;
1158         for ( pIt1 = eData.myPoints.begin(); pIt1 != eData.myPoints.end();  )
1159         {
1160           double distToFace =
1161             triaSearcher->GetMinDistInSphere( pIt1->myP, pIt1->mySegSize, isAdjFace, avoidPnt );
1162           double allowedSize = Max( minSize, distToFace*( 1. + grading ));
1163           if ( 1.1 * allowedSize < pIt1->mySegSize  )
1164           {
1165             sizeDecreased = true;
1166             sizeTree.SetSize( pIt1->myP, allowedSize );
1167             // cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() )
1168             //      << "\t SetSize " << allowedSize << " at "
1169             //      << pIt1->myP.X() <<", "<< pIt1->myP.Y()<<", "<<pIt1->myP.Z() << endl;
1170             pIt2 = pIt1;
1171             if ( --pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
1172               sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1173             pIt2 = pIt1;
1174             if ( ++pIt2 != eData.myPoints.end() && pIt2->mySegSize > allowedSize )
1175               sizeTree.SetSize( eData.myC3d.Value( 0.6*pIt2->myU + 0.4*pIt1->myU ), allowedSize );
1176             pIt1->mySegSize = allowedSize;
1177           }
1178           ++pIt1;
1179           if ( & (*pIt1) == & eData.Last() )
1180             avoidPnt = & eData.Last().myP;
1181           else
1182             avoidPnt = NULL;
1183
1184           if ( iLoop > 20 )
1185           {
1186 #ifdef _DEBUG_
1187             cout << "Infinite loop in AdaptiveAlgo::Compute()" << endl;
1188 #endif
1189             sizeDecreased = false;
1190             break;
1191           }
1192         }
1193       } // while ( sizeDecreased )
1194     } // loop on myEdges
1195
1196     //*theProgress = 0.3 + 0.3 * iF / double( faceMap.Extent() );
1197
1198   } // loop on faceMap
1199
1200   return makeSegments();
1201 }
1202
1203 //================================================================================
1204 /*!
1205  * \brief Create segments
1206  */
1207 //================================================================================
1208
1209 bool AdaptiveAlgo::makeSegments()
1210 {
1211   SMESH_HypoFilter quadHyp( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
1212   _quadraticMesh = myMesh->GetHypothesis( myEdges[0].Edge(), quadHyp, /*andAncestors=*/true );
1213
1214   SMESH_MesherHelper helper( *myMesh );
1215   helper.SetIsQuadratic( _quadraticMesh );
1216
1217   vector< double > nbSegs, params;
1218
1219   for ( int iE = 0; iE < myEdges.size(); ++iE )
1220   {
1221     EdgeData& eData = myEdges[ iE ];
1222
1223     // estimate roughly min segement size on the EDGE
1224     double edgeMinSize = myHyp->GetMaxSize();
1225     EdgeData::TPntIter pIt1 = eData.myPoints.begin();
1226     for ( ; pIt1 != eData.myPoints.end(); ++pIt1 )
1227       edgeMinSize = Min( edgeMinSize, mySizeTree->GetSize( pIt1->myP ));
1228
1229     const double f = eData.myC3d.FirstParameter(), l = eData.myC3d.LastParameter();
1230     const double parLen = l - f;
1231     const int  nbDivSeg = 5;
1232     int           nbDiv = Max( 1, int ( eData.myLength / edgeMinSize * nbDivSeg ));
1233
1234     // compute nb of segments
1235     bool toRecompute = true;
1236     double maxSegSize = 0;
1237     //cout << "E " << theMesh.GetMeshDS()->ShapeToIndex( eData.Edge() ) << endl;
1238     while ( toRecompute ) // recompute if segment size at some point is less than edgeMinSize/nbDivSeg
1239     {
1240       nbSegs.resize( nbDiv + 1 );
1241       nbSegs[0] = 0;
1242       toRecompute = false;
1243
1244       gp_Pnt p1 = eData.First().myP, p2, pDiv = p1;
1245       for ( size_t i = 1, segCount = 1; i < nbSegs.size(); ++i )
1246       {
1247         p2 = eData.myC3d.Value( f + parLen * i / nbDiv );
1248         double locSize = Min( mySizeTree->GetSize( p2 ), myHyp->GetMaxSize() );
1249         double nb      = p1.Distance( p2 ) / locSize;
1250         // if ( nbSegs.size() < 30 )
1251         //   cout << "locSize " << locSize << " nb " << nb << endl;
1252         if ( nb > 1. )
1253         {
1254           toRecompute = true;
1255           edgeMinSize = locSize;
1256           nbDiv = int ( eData.myLength / edgeMinSize * nbDivSeg );
1257           break;
1258         }
1259         nbSegs[i] = nbSegs[i-1] + nb;
1260         p1 = p2;
1261         if ( nbSegs[i] >= segCount )
1262         {
1263           maxSegSize = Max( maxSegSize, pDiv.Distance( p2 ));
1264           pDiv = p2;
1265           ++segCount;
1266         }
1267       }
1268     }
1269
1270     // compute parameters of nodes
1271     int nbSegFinal = Max( 1, int(floor( nbSegs.back() + 0.5 )));
1272     double fact = nbSegFinal / nbSegs.back();
1273     if ( maxSegSize / fact > myHyp->GetMaxSize() )
1274       fact = ++nbSegFinal / nbSegs.back();
1275     //cout << "nbSegs.back() " << nbSegs.back() << " nbSegFinal " << nbSegFinal << endl;
1276     params.clear();
1277     for ( int i = 0, segCount = 1; segCount < nbSegFinal; ++segCount )
1278     {
1279       while ( nbSegs[i] * fact < segCount )
1280         ++i;
1281       if ( i < nbDiv )
1282       {
1283         double d = i - ( nbSegs[i] - segCount/fact ) / ( nbSegs[i] - nbSegs[i-1] );
1284         params.push_back( f + parLen * d / nbDiv );
1285         //params.push_back( f + parLen * i / nbDiv );
1286       }
1287       else
1288         break;
1289     }
1290     // get nodes on VERTEXes
1291     TopoDS_Vertex vf = helper.IthVertex( 0, eData.Edge(), false );
1292     TopoDS_Vertex vl = helper.IthVertex( 1, eData.Edge(), false );
1293     myMesh->GetSubMesh( vf )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1294     myMesh->GetSubMesh( vl )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1295     const SMDS_MeshNode * nf = VertexNode( vf, myMesh->GetMeshDS() );
1296     const SMDS_MeshNode * nl = VertexNode( vl, myMesh->GetMeshDS() );
1297     if ( !nf || !nl )
1298       return error("No node on vertex");
1299
1300     // create segments
1301     helper.SetSubShape( eData.Edge() );
1302     helper.SetElementsOnShape( true );
1303     const int ID = 0;
1304     const SMDS_MeshNode *n1 = nf, *n2;
1305     for ( size_t i = 0; i < params.size(); ++i, n1 = n2 )
1306     {
1307       gp_Pnt p2 = eData.myC3d.Value( params[i] );
1308       n2 = helper.AddNode( p2.X(), p2.Y(), p2.Z(), ID, params[i] );
1309       helper.AddEdge( n1, n2, ID, /*force3d=*/false );
1310     }
1311     helper.AddEdge( n1, nl, ID, /*force3d=*/false );
1312
1313     eData.myPoints.clear();
1314
1315     //*theProgress = 0.6 + 0.4 * iE / double( myEdges.size() );
1316     if ( _computeCanceled )
1317       return false;
1318
1319   } // loop on EDGEs
1320
1321   SMESHUtils::FreeVector( myEdges );
1322
1323   return true;
1324 }
1325
1326 //================================================================================
1327 /*!
1328  * \brief Predict number of segments on all given EDGEs
1329  */
1330 //================================================================================
1331
1332 bool AdaptiveAlgo::Evaluate(SMESH_Mesh &         theMesh,
1333                             const TopoDS_Shape & theShape,
1334                             MapShapeNbElems&     theResMap)
1335 {
1336   // initialize fields of inherited StdMeshers_Regular_1D
1337   StdMeshers_Regular_1D::_hypType = DEFLECTION;
1338   StdMeshers_Regular_1D::_value[ DEFLECTION_IND ] = myHyp->GetDeflection();
1339
1340   TopExp_Explorer edExp( theShape, TopAbs_EDGE );
1341
1342   for ( ; edExp.More(); edExp.Next() )
1343   {
1344     const TopoDS_Edge & edge = TopoDS::Edge( edExp.Current() );
1345     StdMeshers_Regular_1D::Evaluate( theMesh, theShape, theResMap );
1346   }
1347   return true;
1348 }
1349