src/StdMeshers/StdMeshers_Adaptive1D.cxx

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