src/StdMeshers/StdMeshers_Regular_1D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : StdMeshers_Regular_1D.cxx
  25 //           Moved here from SMESH_Regular_1D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //  $Header$
  29
  30 using namespace std;
  31
  32 #include "StdMeshers_Regular_1D.hxx"
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_Mesh.hxx"
  35
  36 #include "StdMeshers_LocalLength.hxx"
  37 #include "StdMeshers_NumberOfSegments.hxx"
  38 #include "StdMeshers_Arithmetic1D.hxx"
  39 #include "StdMeshers_StartEndLength.hxx"
  40 #include "StdMeshers_Deflection1D.hxx"
  41 #include <StdMeshers_AutomaticLength.hxx>
  42
  43 #include "SMDS_MeshElement.hxx"
  44 #include "SMDS_MeshNode.hxx"
  45 #include "SMDS_EdgePosition.hxx"
  46 #include "SMESH_subMesh.hxx"
  47
  48 #include "Utils_SALOME_Exception.hxx"
  49 #include "utilities.h"
  50
  51 #include <BRep_Tool.hxx>
  52 #include <TopoDS_Edge.hxx>
  53 #include <TopoDS_Shape.hxx>
  54 #include <TopTools_ListIteratorOfListOfShape.hxx>
  55 #include <GeomAdaptor_Curve.hxx>
  56 #include <GCPnts_AbscissaPoint.hxx>
  57 #include <GCPnts_UniformAbscissa.hxx>
  58 #include <GCPnts_UniformDeflection.hxx>
  59 #include <Standard_ErrorHandler.hxx>
  60 #include <Precision.hxx>
  61 #include <Expr_GeneralExpression.hxx>
  62 #include <Expr_NamedUnknown.hxx>
  63 #include <Expr_Array1OfNamedUnknown.hxx>
  64 #include <TColStd_Array1OfReal.hxx>
  65 #include <ExprIntrp_GenExp.hxx>
  66
  67 #include <CASCatch_CatchSignals.hxx>
  68 #include <CASCatch_Failure.hxx>
  69 #include <CASCatch_ErrorHandler.hxx>
  70 #include <OSD.hxx>
  71 #include <math_GaussSingleIntegration.hxx>
  72
  73 #include <string>
  74 #include <math.h>
  75
  76 //=============================================================================
  77 /*!
  78  *
  79  */
  80 //=============================================================================
  81
  82 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
  83         SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
  84 {
  85         MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
  86         _name = "Regular_1D";
  87         _shapeType = (1 << TopAbs_EDGE);
  88
  89         _compatibleHypothesis.push_back("LocalLength");
  90         _compatibleHypothesis.push_back("NumberOfSegments");
  91         _compatibleHypothesis.push_back("StartEndLength");
  92         _compatibleHypothesis.push_back("Deflection1D");
  93         _compatibleHypothesis.push_back("Arithmetic1D");
  94         _compatibleHypothesis.push_back("AutomaticLength");
  95 }
  96
  97 //=============================================================================
  98 /*!
  99  *
 100  */
 101 //=============================================================================
 102
 103 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
 104 {
 105 }
 106
 107 //=============================================================================
 108 /*!
 109  *
 110  */
 111 //=============================================================================
 112
 113 bool StdMeshers_Regular_1D::CheckHypothesis
 114                          (SMESH_Mesh& aMesh,
 115                           const TopoDS_Shape& aShape,
 116                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 117 {
 118   _hypType = NONE;
 119
 120   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 121   if (hyps.size() == 0)
 122   {
 123     aStatus = SMESH_Hypothesis::HYP_MISSING;
 124     return false;  // can't work without a hypothesis
 125   }
 126
 127   // use only the first hypothesis
 128   const SMESHDS_Hypothesis *theHyp = hyps.front();
 129
 130   string hypName = theHyp->GetName();
 131
 132   if (hypName == "LocalLength")
 133   {
 134     const StdMeshers_LocalLength * hyp =
 135       dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
 136     ASSERT(hyp);
 137     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength();
 138     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 139     _hypType = LOCAL_LENGTH;
 140     aStatus = SMESH_Hypothesis::HYP_OK;
 141   }
 142
 143   else if (hypName == "NumberOfSegments")
 144   {
 145     const StdMeshers_NumberOfSegments * hyp =
 146       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
 147     ASSERT(hyp);
 148     _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
 149     ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
 150     _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
 151     switch (_ivalue[ DISTR_TYPE_IND ])
 152     {
 153     case StdMeshers_NumberOfSegments::DT_Scale:
 154       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
 155       break;
 156     case StdMeshers_NumberOfSegments::DT_TabFunc:
 157       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
 158       break;
 159     case StdMeshers_NumberOfSegments::DT_ExprFunc:
 160       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
 161       break;
 162     case StdMeshers_NumberOfSegments::DT_Regular:
 163       break;
 164     default:
 165       ASSERT(0);
 166       break;
 167     }
 168     if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
 169         _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
 170       _ivalue[ EXP_MODE_IND ] = (int) hyp->IsExponentMode();
 171     _hypType = NB_SEGMENTS;
 172     aStatus = SMESH_Hypothesis::HYP_OK;
 173   }
 174
 175   else if (hypName == "Arithmetic1D")
 176   {
 177     const StdMeshers_Arithmetic1D * hyp =
 178       dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
 179     ASSERT(hyp);
 180     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 181     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 182     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 183     _hypType = ARITHMETIC_1D;
 184     aStatus = SMESH_Hypothesis::HYP_OK;
 185   }
 186
 187   else if (hypName == "StartEndLength")
 188   {
 189     const StdMeshers_StartEndLength * hyp =
 190       dynamic_cast <const StdMeshers_StartEndLength * >(theHyp);
 191     ASSERT(hyp);
 192     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 193     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 194     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 195     _hypType = BEG_END_LENGTH;
 196     aStatus = SMESH_Hypothesis::HYP_OK;
 197   }
 198
 199   else if (hypName == "Deflection1D")
 200   {
 201     const StdMeshers_Deflection1D * hyp =
 202       dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
 203     ASSERT(hyp);
 204     _value[ DEFLECTION_IND ] = hyp->GetDeflection();
 205     ASSERT( _value[ DEFLECTION_IND ] > 0 );
 206     _hypType = DEFLECTION;
 207     aStatus = SMESH_Hypothesis::HYP_OK;
 208   }
 209
 210   else if (hypName == "AutomaticLength")
 211   {
 212     StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
 213       (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
 214     ASSERT(hyp);
 215     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 216     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 217     _hypType = LOCAL_LENGTH;
 218     aStatus = SMESH_Hypothesis::HYP_OK;
 219   }
 220   else
 221     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 222
 223   return ( _hypType != NONE );
 224 }
 225
 226 //=======================================================================
 227 //function : compensateError
 228 //purpose  : adjust theParams so that the last segment length == an
 229 //=======================================================================
 230
 231 static void compensateError(double a1, double an,
 232                             double U1, double Un,
 233                             double             length,
 234                             GeomAdaptor_Curve& C3d,
 235                             list<double> &     theParams)
 236 {
 237   int i, nPar = theParams.size();
 238   if ( a1 + an < length && nPar > 1 )
 239   {
 240     list<double>::reverse_iterator itU = theParams.rbegin();
 241     double Ul = *itU++;
 242     // dist from the last point to the edge end <Un>, it should be equal <an>
 243     double Ln = GCPnts_AbscissaPoint::Length( C3d, Ul, Un );
 244     double dLn = an - Ln; // error of <an>
 245     if ( Abs( dLn ) <= Precision::Confusion() )
 246       return;
 247     double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
 248     double dUn = dLn * Abs( Un - U1 ) / length; // parametric error of <an>
 249     if ( dUn < 0.5 * dU ) { // last segment is a bit shorter than it should
 250       dUn = -dUn; // move the last parameter to the edge beginning
 251     }
 252     else {  // last segment is much shorter than it should -> remove the last param and
 253       theParams.pop_back(); nPar--; // move the rest points toward the edge end
 254       Ln = GCPnts_AbscissaPoint::Length( C3d, theParams.back(), Un );
 255       dUn = ( an - Ln ) * Abs( Un - U1 ) / length;
 256       if ( dUn < 0.5 * dU )
 257         dUn = -dUn;
 258     }
 259     if ( U1 > Un )
 260       dUn = -dUn;
 261     double q  = dUn / ( nPar - 1 );
 262     for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
 263       (*itU) += dUn;
 264       dUn -= q;
 265     }
 266   }
 267 }
 268
 269 class Function
 270 {
 271 public:
 272   Function( const bool exp )
 273   : myExp( exp )
 274   {
 275   }
 276
 277   virtual ~Function()
 278   {
 279   }
 280
 281   virtual bool   value( const double, double& f )
 282   {
 283     if( myExp )
 284       f = pow( 10, f );
 285     return true;
 286   }
 287   virtual double integral( const double, const double ) = 0;
 288
 289 private:
 290   bool myExp;
 291 };
 292
 293 class FunctionIntegral : public Function
 294 {
 295 public:
 296   FunctionIntegral( Function*, const double );
 297   virtual ~FunctionIntegral();
 298   virtual bool   value( const double, double& );
 299   virtual double integral( const double, const double );
 300
 301 private:
 302   Function* myFunc;
 303   double    myStart;
 304 };
 305
 306 FunctionIntegral::FunctionIntegral( Function* f, const double st )
 307 : Function( false )
 308 {
 309   myFunc = f;
 310   myStart = st;
 311 }
 312
 313 FunctionIntegral::~FunctionIntegral()
 314 {
 315 }
 316
 317 bool FunctionIntegral::value( const double t, double& f )
 318 {
 319   f = myFunc ? myFunc->integral( myStart, t ) : 0;
 320   return myFunc!=0 && Function::value( t, f );
 321 }
 322
 323 double FunctionIntegral::integral( const double, const double )
 324 {
 325   return 0;
 326 }
 327
 328 class FunctionTable : public Function
 329 {
 330 public:
 331   FunctionTable( const std::vector<double>&, const bool );
 332   virtual ~FunctionTable();
 333   virtual bool   value( const double, double& );
 334   virtual double integral( const double, const double );
 335
 336 private:
 337   bool    findBounds( const double, int&, int& ) const;
 338
 339   //integral from x[i] to x[i+1]
 340   double  integral( const int i );
 341
 342   //integral from x[i] to x[i]+d
 343   //warning: function is presented as linear on interaval from x[i] to x[i]+d,
 344   //         for correct result d must be >=0 and <=x[i+1]-x[i]
 345   double  integral( const int i, const double d );
 346
 347 private:
 348   std::vector<double>  myData;
 349 };
 350
 351 FunctionTable::FunctionTable( const std::vector<double>& data, const bool exp )
 352 : Function( exp )
 353 {
 354   myData = data;
 355 }
 356
 357 FunctionTable::~FunctionTable()
 358 {
 359 }
 360
 361 bool FunctionTable::value( const double t, double& f )
 362 {
 363   int i1, i2;
 364   if( !findBounds( t, i1, i2 ) )
 365     return false;
 366
 367   double
 368     x1 = myData[2*i1], y1 = myData[2*i1+1],
 369     x2 = myData[2*i2], y2 = myData[2*i2+1];
 370
 371   Function::value( x1, y1 );
 372   Function::value( x2, y2 );
 373
 374   f = y1 + ( y2-y1 ) * ( t-x1 ) / ( x2-x1 );
 375   return true;
 376 }
 377
 378 double FunctionTable::integral( const int i )
 379 {
 380   if( i>=0 && i<myData.size()-1 )
 381     return integral( i, myData[2*(i+1)]-myData[2*i] );
 382   else
 383     return 0;
 384 }
 385
 386 double FunctionTable::integral( const int i, const double d )
 387 {
 388   double f, res = 0.0;
 389   if( value( myData[2*i]+d, f ) )
 390     res = ( myData[2*i] + f ) / 2.0 * d;
 391
 392   return res;
 393 }
 394
 395 double FunctionTable::integral( const double a, const double b )
 396 {
 397   int x1s, x1f, x2s, x2f;
 398   findBounds( a, x1s, x1f );
 399   findBounds( b, x2s, x2f );
 400   double J = 0;
 401   for( int i=x1s; i<x2s; i++ )
 402     J+=integral( i );
 403   J-=integral( x1s, a-myData[2*x1s] );
 404   J+=integral( x2s, b-myData[2*x2s] );
 405   return J;
 406 }
 407
 408 bool FunctionTable::findBounds( const double x, int& x_ind_1, int& x_ind_2 ) const
 409 {
 410   int n = myData.size();
 411   if( n==0 || x<myData[0] )
 412   {
 413     x_ind_1 = x_ind_2 = 0;
 414     return false;
 415   }
 416
 417   for( int i=0; i<n-1; i++ )
 418     if( myData[2*i]<=x && x<=myData[2*(i+1)] )
 419     {
 420       x_ind_1 = i;
 421       x_ind_2 = i+1;
 422       return true;
 423     }
 424   x_ind_1 = n-1;
 425   x_ind_2 = n-1;
 426   return false;
 427 }
 428
 429
 430
 431 class FunctionExpr : public Function, public math_Function
 432 {
 433 public:
 434   FunctionExpr( const char*, const bool );
 435   virtual ~FunctionExpr();
 436   virtual Standard_Boolean Value( Standard_Real, Standard_Real& );
 437   virtual bool   value( const double, double& );  //inherited from Function
 438   virtual double integral( const double, const double );
 439
 440 private:
 441   Handle(ExprIntrp_GenExp)    myExpr;
 442   Expr_Array1OfNamedUnknown   myVars;
 443   TColStd_Array1OfReal        myValues;
 444 };
 445
 446 FunctionExpr::FunctionExpr( const char* str, const bool exp )
 447 : Function( exp ),
 448   myVars( 1, 1 ),
 449   myValues( 1, 1 )
 450 {
 451   myExpr = ExprIntrp_GenExp::Create();
 452   myExpr->Process( ( Standard_CString )str );
 453   if( !myExpr->IsDone() )
 454     myExpr.Nullify();
 455
 456   myVars.ChangeValue( 1 ) = new Expr_NamedUnknown( "t" );
 457 }
 458
 459 FunctionExpr::~FunctionExpr()
 460 {
 461 }
 462
 463 Standard_Boolean FunctionExpr::Value( Standard_Real T, Standard_Real& F )
 464 {
 465   double f;
 466   Standard_Boolean res = value( T, f );
 467   F = f;
 468   return res;
 469 }
 470
 471 bool FunctionExpr::value( const double t, double& f )
 472 {
 473   if( myExpr.IsNull() )
 474     return false;
 475
 476   CASCatch_CatchSignals aCatchSignals;
 477   aCatchSignals.Activate();
 478
 479   myValues.ChangeValue( 1 ) = t;
 480   bool ok = true;
 481   CASCatch_TRY {
 482     f = myExpr->Expression()->Evaluate( myVars, myValues );
 483   }
 484   CASCatch_CATCH(CASCatch_Failure) {
 485     aCatchSignals.Deactivate();
 486     Handle(CASCatch_Failure) aFail = CASCatch_Failure::Caught();
 487     f = 0.0;
 488   }
 489
 490   aCatchSignals.Deactivate();
 491   ok = Function::value( t, f ) && ok;
 492   return ok;
 493 }
 494
 495 double FunctionExpr::integral( const double a, const double b )
 496 {
 497   double res = 0.0;
 498   CASCatch_TRY
 499   {
 500     math_GaussSingleIntegration _int( *this, a, b, 20 );
 501     if( _int.IsDone() )
 502       res = _int.Value();
 503   }
 504   CASCatch_CATCH(CASCatch_Failure)
 505   {
 506     res = 0.0;
 507     MESSAGE( "Exception in integral calculating" );
 508   }
 509   return res;
 510 }
 511
 512
 513
 514
 515
 516
 517
 518 double dihotomySolve( Function& f, const double val, const double _start, const double _fin, const double eps, bool& ok )
 519 {
 520   double start = _start, fin = _fin, start_val, fin_val; bool ok1, ok2;
 521   ok1 = f.value( start, start_val );
 522   ok2 = f.value( fin, fin_val );
 523
 524   if( !ok1 || !ok2 )
 525   {
 526     ok = false;
 527     return 0.0;
 528   }
 529
 530   bool start_pos = start_val>=val, fin_pos = fin_val>=val;
 531   ok = true;
 532
 533   while( fin-start>eps )
 534   {
 535     double mid = ( start+fin )/2.0, mid_val;
 536     ok = f.value( mid, mid_val );
 537     if( !ok )
 538       return 0.0;
 539
 540 //    char buf[1024];
 541 //    sprintf( buf, "start=%f\nfin=%f\nmid_val=%f\n", float( start ), float( fin ), float( mid_val ) );
 542 //    MESSAGE( buf );
 543
 544     bool mid_pos = mid_val>=val;
 545     if( start_pos!=mid_pos )
 546     {
 547       fin_pos = mid_pos;
 548       fin = mid;
 549     }
 550     else if( fin_pos!=mid_pos )
 551     {
 552       start_pos = mid_pos;
 553       start = mid;
 554     }
 555     else
 556       break;
 557   }
 558   return (start+fin)/2.0;
 559 }
 560
 561 static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
 562                                double length, bool theReverse,
 563                                int nbSeg, Function& func,
 564                                list<double>& theParams)
 565 {
 566   OSD::SetSignal( true );
 567   if( nbSeg<=0 )
 568     return false;
 569
 570   MESSAGE( "computeParamByFunc" );
 571
 572   int nbPnt = 1 + nbSeg;
 573   vector<double> x(nbPnt, 0.);
 574
 575   x[0] = 0.0;
 576   double J = func.integral( 0.0, 1.0 ) / nbSeg;
 577   bool ok;
 578   for( int i=1; i<nbSeg; i++ )
 579   {
 580     FunctionIntegral f_int( &func, x[i-1] );
 581     x[i] = dihotomySolve( f_int, J, x[i-1], 1.0, 1E-4, ok );
 582     if( !ok )
 583       return false;
 584   }
 585
 586   x[nbSeg] = 1.0;
 587   MESSAGE( "Points:\n" );
 588   char buf[1024];
 589   for( int i=0; i<=nbSeg; i++ )
 590   {
 591     sprintf(  buf, "%f\n", float(x[i] ) );
 592     MESSAGE( buf );
 593   }
 594
 595
 596
 597   // apply parameters in range [0,1] to the space of the curve
 598   double prevU = first;
 599   double sign = 1.;
 600   if (theReverse)
 601   {
 602     prevU = last;
 603     sign = -1.;
 604   }
 605   for( int i = 1; i < nbSeg; i++ )
 606   {
 607     double curvLength = length * (x[i] - x[i-1]) * sign;
 608     GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
 609     if ( !Discret.IsDone() )
 610       return false;
 611     double U = Discret.Parameter();
 612     if ( U > first && U < last )
 613       theParams.push_back( U );
 614     else
 615       return false;
 616     prevU = U;
 617   }
 618   return true;
 619 }
 620
 621 //=============================================================================
 622 /*!
 623  *
 624  */
 625 //=============================================================================
 626 bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge,
 627                                                       list<double> &     theParams,
 628                                                       const bool         theReverse) const
 629 {
 630   theParams.clear();
 631
 632   double f, l;
 633   Handle(Geom_Curve) Curve = BRep_Tool::Curve(theEdge, f, l);
 634   GeomAdaptor_Curve C3d(Curve);
 635
 636   double length = EdgeLength(theEdge);
 637
 638   switch( _hypType )
 639   {
 640   case LOCAL_LENGTH:
 641   case NB_SEGMENTS: {
 642
 643     double eltSize = 1;
 644     if ( _hypType == LOCAL_LENGTH )
 645     {
 646       // Local Length hypothesis
 647       double nbseg = ceil(length / _value[ BEG_LENGTH_IND ]); // integer sup
 648       if (nbseg <= 0)
 649         nbseg = 1;                        // degenerated edge
 650       eltSize = length / nbseg;
 651     }
 652     else
 653     {
 654       // Number Of Segments hypothesis
 655       switch (_ivalue[ DISTR_TYPE_IND ])
 656       {
 657       case StdMeshers_NumberOfSegments::DT_Scale:
 658         {
 659           double scale = _value[ SCALE_FACTOR_IND ];
 660           if ( theReverse )
 661             scale = 1. / scale;
 662           double alpha = pow( scale , 1.0 / (_ivalue[ NB_SEGMENTS_IND ] - 1));
 663           double factor = (l - f) / (1 - pow( alpha,_ivalue[ NB_SEGMENTS_IND ]));
 664
 665           int i, NbPoints = 1 + _ivalue[ NB_SEGMENTS_IND ];
 666           for ( i = 2; i < NbPoints; i++ )
 667           {
 668             double param = f + factor * (1 - pow(alpha, i - 1));
 669             theParams.push_back( param );
 670           }
 671           return true;
 672         }
 673         break;
 674       case StdMeshers_NumberOfSegments::DT_TabFunc:
 675         {
 676           FunctionTable func(_vvalue[ TAB_FUNC_IND ], (bool)_ivalue[ EXP_MODE_IND ]);
 677           return computeParamByFunc(C3d, f, l, length, theReverse,
 678                                     _ivalue[ NB_SEGMENTS_IND ], func,
 679                                     theParams);
 680         }
 681         break;
 682       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 683         {
 684           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), (bool)_ivalue[ EXP_MODE_IND ]);
 685           return computeParamByFunc(C3d, f, l, length, theReverse,
 686                                     _ivalue[ NB_SEGMENTS_IND ], func,
 687                                     theParams);
 688         }
 689         break;
 690       case StdMeshers_NumberOfSegments::DT_Regular:
 691         eltSize = length / _ivalue[ NB_SEGMENTS_IND ];
 692         break;
 693       default:
 694         return false;
 695       }
 696     }
 697
 698     GCPnts_UniformAbscissa Discret(C3d, eltSize, f, l);
 699     if ( !Discret.IsDone() )
 700       return false;
 701
 702     int NbPoints = Discret.NbPoints();
 703     for ( int i = 2; i < NbPoints; i++ )
 704     {
 705       double param = Discret.Parameter(i);
 706       theParams.push_back( param );
 707     }
 708     return true;
 709   }
 710
 711   case BEG_END_LENGTH: {
 712
 713     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = length
 714
 715     double a1 = _value[ BEG_LENGTH_IND ];
 716     double an = _value[ END_LENGTH_IND ];
 717     double q  = ( length - a1 ) / ( length - an );
 718
 719     double U1 = theReverse ? l : f;
 720     double Un = theReverse ? f : l;
 721     double param = U1;
 722     double eltSize = theReverse ? -a1 : a1;
 723     while ( 1 ) {
 724       // computes a point on a curve <C3d> at the distance <eltSize>
 725       // from the point of parameter <param>.
 726       GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
 727       if ( !Discret.IsDone() ) break;
 728       param = Discret.Parameter();
 729       if ( param > f && param < l )
 730         theParams.push_back( param );
 731       else
 732         break;
 733       eltSize *= q;
 734     }
 735     compensateError( a1, an, U1, Un, length, C3d, theParams );
 736     return true;
 737   }
 738
 739   case ARITHMETIC_1D: {
 740
 741     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = length
 742
 743     double a1 = _value[ BEG_LENGTH_IND ];
 744     double an = _value[ END_LENGTH_IND ];
 745
 746     double  q = ( an - a1 ) / ( 2 *length/( a1 + an ) - 1 );
 747     int     n = int( 1 + ( an - a1 ) / q );
 748
 749     double U1 = theReverse ? l : f;
 750     double Un = theReverse ? f : l;
 751     double param = U1;
 752     double eltSize = a1;
 753     if ( theReverse ) {
 754       eltSize = -eltSize;
 755       q = -q;
 756     }
 757     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 758       // computes a point on a curve <C3d> at the distance <eltSize>
 759       // from the point of parameter <param>.
 760       GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
 761       if ( !Discret.IsDone() ) break;
 762       param = Discret.Parameter();
 763       if ( param > f && param < l )
 764         theParams.push_back( param );
 765       else
 766         break;
 767       eltSize += q;
 768     }
 769     compensateError( a1, an, U1, Un, length, C3d, theParams );
 770
 771     return true;
 772   }
 773
 774   case DEFLECTION: {
 775
 776     GCPnts_UniformDeflection Discret(C3d, _value[ DEFLECTION_IND ], true);
 777     if ( !Discret.IsDone() )
 778       return false;
 779
 780     int NbPoints = Discret.NbPoints();
 781     for ( int i = 2; i < NbPoints; i++ )
 782     {
 783       double param = Discret.Parameter(i);
 784       theParams.push_back( param );
 785     }
 786     return true;
 787
 788   }
 789
 790   default:;
 791   }
 792
 793   return false;
 794 }
 795
 796 //=============================================================================
 797 /*!
 798  *
 799  */
 800 //=============================================================================
 801
 802 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
 803 {
 804   MESSAGE("StdMeshers_Regular_1D::Compute");
 805
 806   if ( _hypType == NONE )
 807     return false;
 808
 809   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 810   aMesh.GetSubMesh(aShape);
 811
 812   const TopoDS_Edge & EE = TopoDS::Edge(aShape);
 813   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
 814   int shapeID = meshDS->ShapeToIndex( E );
 815
 816   double f, l;
 817   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
 818
 819   TopoDS_Vertex VFirst, VLast;
 820   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
 821
 822   ASSERT(!VFirst.IsNull());
 823   SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
 824   if (!lid->more())
 825   {
 826     MESSAGE (" NO NODE BUILT ON VERTEX ");
 827     return false;
 828   }
 829   const SMDS_MeshNode * idFirst = lid->next();
 830
 831   ASSERT(!VLast.IsNull());
 832   lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
 833   if (!lid->more())
 834   {
 835     MESSAGE (" NO NODE BUILT ON VERTEX ");
 836     return false;
 837   }
 838   const SMDS_MeshNode * idLast = lid->next();
 839
 840   if (!Curve.IsNull())
 841   {
 842     list< double > params;
 843     bool reversed = false;
 844     if ( !_mainEdge.IsNull() )
 845       reversed = aMesh.IsReversedInChain( EE, _mainEdge );
 846     try {
 847       if ( ! computeInternalParameters( E, params, reversed ))
 848         return false;
 849     }
 850     catch ( Standard_Failure ) {
 851       return false;
 852     }
 853
 854     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
 855     // only internal nodes receive an edge position with param on curve
 856
 857     const SMDS_MeshNode * idPrev = idFirst;
 858
 859     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 860     {
 861       double param = *itU;
 862       gp_Pnt P = Curve->Value(param);
 863
 864       //Add the Node in the DataStructure
 865       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 866       meshDS->SetNodeOnEdge(node, shapeID, param);
 867
 868       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 869       meshDS->SetMeshElementOnShape(edge, shapeID);
 870       idPrev = node;
 871     }
 872     SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
 873     meshDS->SetMeshElementOnShape(edge, shapeID);
 874   }
 875   else
 876   {
 877     // Edge is a degenerated Edge : We put n = 5 points on the edge.
 878     int NbPoints = 5;
 879     BRep_Tool::Range(E, f, l);
 880     double du = (l - f) / (NbPoints - 1);
 881     //MESSAGE("************* Degenerated edge! *****************");
 882
 883     TopoDS_Vertex V1, V2;
 884     TopExp::Vertices(E, V1, V2);
 885     gp_Pnt P = BRep_Tool::Pnt(V1);
 886
 887     const SMDS_MeshNode * idPrev = idFirst;
 888     for (int i = 2; i < NbPoints; i++)
 889     {
 890       double param = f + (i - 1) * du;
 891       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 892       meshDS->SetNodeOnEdge(node, shapeID, param);
 893
 894       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 895       meshDS->SetMeshElementOnShape(edge, shapeID);
 896       idPrev = node;
 897     }
 898     SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
 899     meshDS->SetMeshElementOnShape(edge, shapeID);
 900   }
 901   return true;
 902 }
 903
 904 //=============================================================================
 905 /*!
 906  *  See comments in SMESH_Algo.cxx
 907  */
 908 //=============================================================================
 909
 910 const list <const SMESHDS_Hypothesis *> & StdMeshers_Regular_1D::GetUsedHypothesis(
 911         SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
 912 {
 913   _usedHypList.clear();
 914   _usedHypList = GetAppliedHypothesis(aMesh, aShape);   // copy
 915   int nbHyp = _usedHypList.size();
 916   _mainEdge.Nullify();
 917   if (nbHyp == 0)
 918   {
 919     // Check, if propagated from some other edge
 920     if (aShape.ShapeType() == TopAbs_EDGE &&
 921         aMesh.IsPropagatedHypothesis(aShape, _mainEdge))
 922     {
 923       // Propagation of 1D hypothesis from <aMainEdge> on this edge
 924       //_usedHypList = GetAppliedHypothesis(aMesh, _mainEdge);  // copy
 925       // use a general method in order not to nullify _mainEdge
 926       _usedHypList = SMESH_Algo::GetUsedHypothesis(aMesh, _mainEdge);   // copy
 927       nbHyp = _usedHypList.size();
 928     }
 929   }
 930   if (nbHyp == 0)
 931   {
 932     TopTools_ListIteratorOfListOfShape ancIt( aMesh.GetAncestors( aShape ));
 933     for (; ancIt.More(); ancIt.Next())
 934     {
 935       const TopoDS_Shape& ancestor = ancIt.Value();
 936       _usedHypList = GetAppliedHypothesis(aMesh, ancestor);     // copy
 937       nbHyp = _usedHypList.size();
 938       if (nbHyp == 1)
 939         break;
 940     }
 941   }
 942   if (nbHyp > 1)
 943     _usedHypList.clear();       //only one compatible hypothesis allowed
 944   return _usedHypList;
 945 }
 946
 947 //=============================================================================
 948 /*!
 949  *
 950  */
 951 //=============================================================================
 952
 953 ostream & StdMeshers_Regular_1D::SaveTo(ostream & save)
 954 {
 955   return save;
 956 }
 957
 958 //=============================================================================
 959 /*!
 960  *
 961  */
 962 //=============================================================================
 963
 964 istream & StdMeshers_Regular_1D::LoadFrom(istream & load)
 965 {
 966   return load;
 967 }
 968
 969 //=============================================================================
 970 /*!
 971  *
 972  */
 973 //=============================================================================
 974
 975 ostream & operator <<(ostream & save, StdMeshers_Regular_1D & hyp)
 976 {
 977   return hyp.SaveTo( save );
 978 }
 979
 980 //=============================================================================
 981 /*!
 982  *
 983  */
 984 //=============================================================================
 985
 986 istream & operator >>(istream & load, StdMeshers_Regular_1D & hyp)
 987 {
 988   return hyp.LoadFrom( load );
 989 }