src/StdMeshers/StdMeshers_Adaptive1D.cxx

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