src/StdMeshers/StdMeshers_Adaptive1D.cxx

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