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