src/StdMeshers/StdMeshers_Adaptive1D.cxx

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