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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  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 #include "StdMeshers_Regular_1D.hxx"
  31 #include "StdMeshers_Distribution.hxx"
  32
  33 #include "StdMeshers_LocalLength.hxx"
  34 #include "StdMeshers_NumberOfSegments.hxx"
  35 #include "StdMeshers_Arithmetic1D.hxx"
  36 #include "StdMeshers_StartEndLength.hxx"
  37 #include "StdMeshers_Deflection1D.hxx"
  38 #include "StdMeshers_AutomaticLength.hxx"
  39 #include "StdMeshers_SegmentLengthAroundVertex.hxx"
  40
  41 #include "SMESH_Gen.hxx"
  42 #include "SMESH_Mesh.hxx"
  43 #include "SMESH_HypoFilter.hxx"
  44 #include "SMESH_subMesh.hxx"
  45 #include "SMESH_subMeshEventListener.hxx"
  46 #include "SMESH_Comment.hxx"
  47
  48 #include "SMDS_MeshElement.hxx"
  49 #include "SMDS_MeshNode.hxx"
  50
  51 #include "Utils_SALOME_Exception.hxx"
  52 #include "utilities.h"
  53
  54 #include <BRepAdaptor_Curve.hxx>
  55 #include <BRep_Tool.hxx>
  56 #include <TopoDS_Edge.hxx>
  57 #include <TopExp_Explorer.hxx>
  58 #include <GCPnts_AbscissaPoint.hxx>
  59 #include <GCPnts_UniformAbscissa.hxx>
  60 #include <GCPnts_UniformDeflection.hxx>
  61 #include <Precision.hxx>
  62
  63 #include <string>
  64
  65 using namespace std;
  66
  67 //=============================================================================
  68 /*!
  69  *
  70  */
  71 //=============================================================================
  72
  73 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
  74         SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
  75 {
  76         MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
  77         _name = "Regular_1D";
  78         _shapeType = (1 << TopAbs_EDGE);
  79
  80         _compatibleHypothesis.push_back("LocalLength");
  81         _compatibleHypothesis.push_back("NumberOfSegments");
  82         _compatibleHypothesis.push_back("StartEndLength");
  83         _compatibleHypothesis.push_back("Deflection1D");
  84         _compatibleHypothesis.push_back("Arithmetic1D");
  85         _compatibleHypothesis.push_back("AutomaticLength");
  86
  87         _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
  88 }
  89
  90 //=============================================================================
  91 /*!
  92  *
  93  */
  94 //=============================================================================
  95
  96 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
  97 {
  98 }
  99
 100 //=============================================================================
 101 /*!
 102  *
 103  */
 104 //=============================================================================
 105
 106 bool StdMeshers_Regular_1D::CheckHypothesis
 107                          (SMESH_Mesh&                          aMesh,
 108                           const TopoDS_Shape&                  aShape,
 109                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 110 {
 111   _hypType = NONE;
 112   _quadraticMesh = false;
 113
 114   const bool ignoreAuxiliaryHyps = false;
 115   const list <const SMESHDS_Hypothesis * > & hyps =
 116     GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
 117
 118   // find non-auxiliary hypothesis
 119   const SMESHDS_Hypothesis *theHyp = 0;
 120   list <const SMESHDS_Hypothesis * >::const_iterator h = hyps.begin();
 121   for ( ; h != hyps.end(); ++h ) {
 122     if ( static_cast<const SMESH_Hypothesis*>(*h)->IsAuxiliary() ) {
 123       if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 )
 124         _quadraticMesh = true;
 125     }
 126     else {
 127       if ( !theHyp )
 128         theHyp = *h; // use only the first non-auxiliary hypothesis
 129     }
 130   }
 131
 132   if ( !theHyp )
 133   {
 134     aStatus = SMESH_Hypothesis::HYP_MISSING;
 135     return false;  // can't work without a hypothesis
 136   }
 137
 138   string hypName = theHyp->GetName();
 139
 140   if (hypName == "LocalLength")
 141   {
 142     const StdMeshers_LocalLength * hyp =
 143       dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
 144     ASSERT(hyp);
 145     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength();
 146     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 147     _hypType = LOCAL_LENGTH;
 148     aStatus = SMESH_Hypothesis::HYP_OK;
 149   }
 150
 151   else if (hypName == "NumberOfSegments")
 152   {
 153     const StdMeshers_NumberOfSegments * hyp =
 154       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
 155     ASSERT(hyp);
 156     _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
 157     ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
 158     _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
 159     switch (_ivalue[ DISTR_TYPE_IND ])
 160     {
 161     case StdMeshers_NumberOfSegments::DT_Scale:
 162       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
 163       break;
 164     case StdMeshers_NumberOfSegments::DT_TabFunc:
 165       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
 166       break;
 167     case StdMeshers_NumberOfSegments::DT_ExprFunc:
 168       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
 169       break;
 170     case StdMeshers_NumberOfSegments::DT_Regular:
 171       break;
 172     default:
 173       ASSERT(0);
 174       break;
 175     }
 176     if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
 177         _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
 178         _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
 179     _hypType = NB_SEGMENTS;
 180     aStatus = SMESH_Hypothesis::HYP_OK;
 181   }
 182
 183   else if (hypName == "Arithmetic1D")
 184   {
 185     const StdMeshers_Arithmetic1D * hyp =
 186       dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
 187     ASSERT(hyp);
 188     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 189     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 190     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 191     _hypType = ARITHMETIC_1D;
 192     aStatus = SMESH_Hypothesis::HYP_OK;
 193   }
 194
 195   else if (hypName == "StartEndLength")
 196   {
 197     const StdMeshers_StartEndLength * hyp =
 198       dynamic_cast <const StdMeshers_StartEndLength * >(theHyp);
 199     ASSERT(hyp);
 200     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 201     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 202     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 203     _hypType = BEG_END_LENGTH;
 204     aStatus = SMESH_Hypothesis::HYP_OK;
 205   }
 206
 207   else if (hypName == "Deflection1D")
 208   {
 209     const StdMeshers_Deflection1D * hyp =
 210       dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
 211     ASSERT(hyp);
 212     _value[ DEFLECTION_IND ] = hyp->GetDeflection();
 213     ASSERT( _value[ DEFLECTION_IND ] > 0 );
 214     _hypType = DEFLECTION;
 215     aStatus = SMESH_Hypothesis::HYP_OK;
 216   }
 217
 218   else if (hypName == "AutomaticLength")
 219   {
 220     StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
 221       (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
 222     ASSERT(hyp);
 223     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 224     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 225     _hypType = LOCAL_LENGTH;
 226     aStatus = SMESH_Hypothesis::HYP_OK;
 227   }
 228   else
 229     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 230
 231   return ( _hypType != NONE );
 232 }
 233
 234 static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
 235                                double length, bool theReverse,
 236                                int nbSeg, Function& func,
 237                                list<double>& theParams)
 238 {
 239   // never do this way
 240   //OSD::SetSignal( true );
 241
 242   if( nbSeg<=0 )
 243     return false;
 244
 245   MESSAGE( "computeParamByFunc" );
 246
 247   int nbPnt = 1 + nbSeg;
 248   vector<double> x(nbPnt, 0.);
 249
 250   if( !buildDistribution( func, 0.0, 1.0, nbSeg, x, 1E-4 ) )
 251      return false;
 252
 253   MESSAGE( "Points:\n" );
 254   char buf[1024];
 255   for( int i=0; i<=nbSeg; i++ )
 256   {
 257     sprintf(  buf, "%f\n", float(x[i] ) );
 258     MESSAGE( buf );
 259   }
 260
 261
 262
 263   // apply parameters in range [0,1] to the space of the curve
 264   double prevU = first;
 265   double sign = 1.;
 266   if (theReverse)
 267   {
 268     prevU = last;
 269     sign = -1.;
 270   }
 271   for( int i = 1; i < nbSeg; i++ )
 272   {
 273     double curvLength = length * (x[i] - x[i-1]) * sign;
 274     GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
 275     if ( !Discret.IsDone() )
 276       return false;
 277     double U = Discret.Parameter();
 278     if ( U > first && U < last )
 279       theParams.push_back( U );
 280     else
 281       return false;
 282     prevU = U;
 283   }
 284   return true;
 285 }
 286
 287
 288 //================================================================================
 289 /*!
 290  * \brief adjust internal node parameters so that the last segment length == an
 291   * \param a1 - the first segment length
 292   * \param an - the last segment length
 293   * \param U1 - the first edge parameter
 294   * \param Un - the last edge parameter
 295   * \param length - the edge length
 296   * \param C3d - the edge curve
 297   * \param theParams - internal node parameters to adjust
 298   * \param adjustNeighbors2an - to adjust length of segments next to the last one
 299   *  and not to remove parameters
 300  */
 301 //================================================================================
 302
 303 static void compensateError(double a1, double an,
 304                             double U1, double Un,
 305                             double            length,
 306                             Adaptor3d_Curve&  C3d,
 307                             list<double> &    theParams,
 308                             bool              adjustNeighbors2an = false)
 309 {
 310   int i, nPar = theParams.size();
 311   if ( a1 + an < length && nPar > 1 )
 312   {
 313     bool reverse = ( U1 > Un );
 314     GCPnts_AbscissaPoint Discret(C3d, reverse ? an : -an, Un);
 315     if ( !Discret.IsDone() )
 316       return;
 317     double Utgt = Discret.Parameter(); // target value of the last parameter
 318     list<double>::reverse_iterator itU = theParams.rbegin();
 319     double Ul = *itU++; // real value of the last parameter
 320     double dUn = Utgt - Ul; // parametric error of <an>
 321     if ( Abs(dUn) <= Precision::Confusion() )
 322       return;
 323     double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
 324     if ( adjustNeighbors2an || Abs(dUn) < 0.5 * dU ) { // last segment is a bit shorter than it should
 325       // move the last parameter to the edge beginning
 326     }
 327     else {  // last segment is much shorter than it should -> remove the last param and
 328       theParams.pop_back(); nPar--; // move the rest points toward the edge end
 329       dUn = Utgt - theParams.back();
 330     }
 331
 332     double q  = dUn / ( nPar - 1 );
 333     if ( !adjustNeighbors2an ) {
 334       for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
 335         (*itU) += dUn;
 336         dUn -= q;
 337       }
 338     }
 339     else {
 340       theParams.back() += dUn;
 341       double sign = reverse ? -1 : 1;
 342       double prevU = theParams.back();
 343       itU = theParams.rbegin();
 344       for ( ++itU, i = 2; i < nPar; ++itU, i++ ) {
 345         double newU = *itU + dUn;
 346         if ( newU*sign < prevU*sign ) {
 347           prevU = *itU = newU;
 348           dUn -= q;
 349         }
 350         else { // set U between prevU and next valid param
 351           list<double>::reverse_iterator itU2 = itU;
 352           ++itU2;
 353           int nb = 2;
 354           while ( (*itU2)*sign > prevU*sign ) {
 355             ++itU2; ++nb;
 356           }
 357           dU = ( *itU2 - prevU ) / nb;
 358           while ( itU != itU2 ) {
 359             *itU += dU; ++itU;
 360           }
 361           break;
 362         }
 363       }
 364     }
 365   }
 366 }
 367
 368 //================================================================================
 369 /*!
 370  * \brief Class used to clean mesh on edges when 0D hyp modified.
 371  * Common approach doesn't work when 0D algo is missing because the 0D hyp is
 372  * considered as not participating in computation whereas it is used by 1D algo.
 373  */
 374 //================================================================================
 375
 376 // struct VertexEventListener : public SMESH_subMeshEventListener
 377 // {
 378 //   VertexEventListener():SMESH_subMeshEventListener(0) // won't be deleted by submesh
 379 //   {}
 380 //   /*!
 381 //    * \brief Clean mesh on edges
 382 //    * \param event - algo_event or compute_event itself (of SMESH_subMesh)
 383 //    * \param eventType - ALGO_EVENT or COMPUTE_EVENT (of SMESH_subMesh)
 384 //    * \param subMesh - the submesh where the event occures
 385 //    */
 386 //   void ProcessEvent(const int event, const int eventType, SMESH_subMesh* subMesh,
 387 //                     EventListenerData*, const SMESH_Hypothesis*)
 388 //   {
 389 //     if ( eventType == SMESH_subMesh::ALGO_EVENT) // all algo events
 390 //     {
 391 //       subMesh->ComputeStateEngine( SMESH_subMesh::MODIF_ALGO_STATE );
 392 //     }
 393 //   }
 394 // }; // struct VertexEventListener
 395
 396 //=============================================================================
 397 /*!
 398  * \brief Sets event listener to vertex submeshes
 399  * \param subMesh - submesh where algo is set
 400  *
 401  * This method is called when a submesh gets HYP_OK algo_state.
 402  * After being set, event listener is notified on each event of a submesh.
 403  */
 404 //=============================================================================
 405
 406 void StdMeshers_Regular_1D::SetEventListener(SMESH_subMesh* subMesh)
 407 {
 408 //   static VertexEventListener listener;
 409 //   SMESH_subMeshIteratorPtr smIt = subMesh->getDependsOnIterator(false,false);
 410 //   while (smIt->more()) {
 411 //     subMesh->SetEventListener( &listener, 0, smIt->next() );
 412 //   }
 413 }
 414
 415 //=============================================================================
 416 /*!
 417  * \brief Do nothing
 418  * \param subMesh - restored submesh
 419  *
 420  * This method is called only if a submesh has HYP_OK algo_state.
 421  */
 422 //=============================================================================
 423
 424 void StdMeshers_Regular_1D::SubmeshRestored(SMESH_subMesh* subMesh)
 425 {
 426 }
 427
 428 //=============================================================================
 429 /*!
 430  * \brief Return StdMeshers_SegmentLengthAroundVertex assigned to vertex
 431  */
 432 //=============================================================================
 433
 434 const StdMeshers_SegmentLengthAroundVertex*
 435 StdMeshers_Regular_1D::getVertexHyp(SMESH_Mesh &          theMesh,
 436                                     const TopoDS_Vertex & theV)
 437 {
 438   static SMESH_HypoFilter filter( SMESH_HypoFilter::HasName("SegmentAroundVertex_0D"));
 439   if ( const SMESH_Hypothesis * h = theMesh.GetHypothesis( theV, filter, true ))
 440   {
 441     SMESH_Algo* algo = const_cast< SMESH_Algo* >( static_cast< const SMESH_Algo* > ( h ));
 442     const list <const SMESHDS_Hypothesis *> & hypList = algo->GetUsedHypothesis( theMesh, theV, 0 );
 443     if ( !hypList.empty() && string("SegmentLengthAroundVertex") == hypList.front()->GetName() )
 444       return static_cast<const StdMeshers_SegmentLengthAroundVertex*>( hypList.front() );
 445   }
 446   return 0;
 447 }
 448
 449 //================================================================================
 450 /*!
 451  * \brief Tune parameters to fit "SegmentLengthAroundVertex" hypothesis
 452   * \param theC3d - wire curve
 453   * \param theLength - curve length
 454   * \param theParameters - internal nodes parameters to modify
 455   * \param theVf - 1st vertex
 456   * \param theVl - 2nd vertex
 457  */
 458 //================================================================================
 459
 460 void StdMeshers_Regular_1D::redistributeNearVertices (SMESH_Mesh &          theMesh,
 461                                                       Adaptor3d_Curve &     theC3d,
 462                                                       double                theLength,
 463                                                       std::list< double > & theParameters,
 464                                                       const TopoDS_Vertex & theVf,
 465                                                       const TopoDS_Vertex & theVl)
 466 {
 467   double f = theC3d.FirstParameter(), l = theC3d.LastParameter();
 468   int nPar = theParameters.size();
 469   for ( int isEnd1 = 0; isEnd1 < 2; ++isEnd1 )
 470   {
 471     const TopoDS_Vertex & V = isEnd1 ? theVf : theVl;
 472     const StdMeshers_SegmentLengthAroundVertex* hyp = getVertexHyp (theMesh, V );
 473     if ( hyp ) {
 474       double vertexLength = hyp->GetLength();
 475       if ( vertexLength > theLength / 2.0 )
 476         continue;
 477       if ( isEnd1 ) { // to have a segment of interest at end of theParameters
 478         theParameters.reverse();
 479         std::swap( f, l );
 480       }
 481       if ( _hypType == NB_SEGMENTS )
 482       {
 483         compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 484       }
 485       else if ( nPar <= 3 )
 486       {
 487         if ( !isEnd1 )
 488           vertexLength = -vertexLength;
 489         GCPnts_AbscissaPoint Discret(theC3d, vertexLength, l);
 490         if ( Discret.IsDone() ) {
 491           if ( nPar == 0 )
 492             theParameters.push_back( Discret.Parameter());
 493           else {
 494             double L = GCPnts_AbscissaPoint::Length( theC3d, theParameters.back(), l);
 495             if ( vertexLength < L / 2.0 )
 496               theParameters.push_back( Discret.Parameter());
 497             else
 498               compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 499           }
 500         }
 501       }
 502       else
 503       {
 504         // recompute params between the last segment and a middle one.
 505         // find size of a middle segment
 506         int nHalf = ( nPar-1 ) / 2;
 507         list< double >::reverse_iterator itU = theParameters.rbegin();
 508         std::advance( itU, nHalf );
 509         double Um = *itU++;
 510         double Lm = GCPnts_AbscissaPoint::Length( theC3d, Um, *itU);
 511         double L = GCPnts_AbscissaPoint::Length( theC3d, *itU, l);
 512         StdMeshers_Regular_1D algo( *this );
 513         algo._hypType = BEG_END_LENGTH;
 514         algo._value[ BEG_LENGTH_IND ] = Lm;
 515         algo._value[ END_LENGTH_IND ] = vertexLength;
 516         double from = *itU, to = l;
 517         if ( isEnd1 ) {
 518           std::swap( from, to );
 519           std::swap( algo._value[ BEG_LENGTH_IND ], algo._value[ END_LENGTH_IND ]);
 520         }
 521         list<double> params;
 522         if ( algo.computeInternalParameters( theC3d, L, from, to, params, false ))
 523         {
 524           if ( isEnd1 ) params.reverse();
 525           while ( 1 + nHalf-- )
 526             theParameters.pop_back();
 527           theParameters.splice( theParameters.end(), params );
 528         }
 529         else
 530         {
 531           compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 532         }
 533       }
 534       if ( isEnd1 )
 535         theParameters.reverse();
 536     }
 537   }
 538 }
 539
 540 //=============================================================================
 541 /*!
 542  *
 543  */
 544 //=============================================================================
 545 bool StdMeshers_Regular_1D::computeInternalParameters(Adaptor3d_Curve& theC3d,
 546                                                       double           theLength,
 547                                                       double           theFirstU,
 548                                                       double           theLastU,
 549                                                       list<double> &   theParams,
 550                                                       const bool       theReverse)
 551 {
 552   theParams.clear();
 553
 554   double f = theFirstU, l = theLastU;
 555
 556   switch( _hypType )
 557   {
 558   case LOCAL_LENGTH:
 559   case NB_SEGMENTS: {
 560
 561     double eltSize = 1;
 562     if ( _hypType == LOCAL_LENGTH )
 563     {
 564       // Local Length hypothesis
 565       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 566       if (nbseg <= 0)
 567         nbseg = 1;                        // degenerated edge
 568       eltSize = theLength / nbseg;
 569     }
 570     else
 571     {
 572       // Number Of Segments hypothesis
 573       int NbSegm = _ivalue[ NB_SEGMENTS_IND ];
 574       if ( NbSegm < 1 )  return false;
 575       if ( NbSegm == 1 ) return true;
 576
 577       switch (_ivalue[ DISTR_TYPE_IND ])
 578       {
 579       case StdMeshers_NumberOfSegments::DT_Scale:
 580         {
 581           double scale = _value[ SCALE_FACTOR_IND ];
 582
 583           if (fabs(scale - 1.0) < Precision::Confusion()) {
 584             // special case to avoid division on zero
 585             for (int i = 1; i < NbSegm; i++) {
 586               double param = f + (l - f) * i / NbSegm;
 587               theParams.push_back( param );
 588             }
 589           } else {
 590             // general case of scale distribution
 591             if ( theReverse )
 592               scale = 1.0 / scale;
 593
 594             double alpha = pow(scale, 1.0 / (NbSegm - 1));
 595             double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
 596
 597             for (int i = 1; i < NbSegm; i++) {
 598               double param = f + factor * (1.0 - pow(alpha, i));
 599               theParams.push_back( param );
 600             }
 601           }
 602           return true;
 603         }
 604         break;
 605       case StdMeshers_NumberOfSegments::DT_TabFunc:
 606         {
 607           FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
 608           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 609                                     _ivalue[ NB_SEGMENTS_IND ], func,
 610                                     theParams);
 611         }
 612         break;
 613       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 614         {
 615           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
 616           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 617                                     _ivalue[ NB_SEGMENTS_IND ], func,
 618                                     theParams);
 619         }
 620         break;
 621       case StdMeshers_NumberOfSegments::DT_Regular:
 622         eltSize = theLength / _ivalue[ NB_SEGMENTS_IND ];
 623         break;
 624       default:
 625         return false;
 626       }
 627     }
 628     GCPnts_UniformAbscissa Discret(theC3d, eltSize, f, l);
 629     if ( !Discret.IsDone() )
 630       return error( dfltErr(), "GCPnts_UniformAbscissa failed");
 631
 632     int NbPoints = Discret.NbPoints();
 633     for ( int i = 2; i < NbPoints; i++ )
 634     {
 635       double param = Discret.Parameter(i);
 636       theParams.push_back( param );
 637     }
 638     compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams ); // for PAL9899
 639     return true;
 640   }
 641
 642   case BEG_END_LENGTH: {
 643
 644     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = theLength
 645
 646     double a1 = _value[ BEG_LENGTH_IND ];
 647     double an = _value[ END_LENGTH_IND ];
 648     double q  = ( theLength - a1 ) / ( theLength - an );
 649
 650     double U1 = theReverse ? l : f;
 651     double Un = theReverse ? f : l;
 652     double param = U1;
 653     double eltSize = theReverse ? -a1 : a1;
 654     while ( 1 ) {
 655       // computes a point on a curve <theC3d> at the distance <eltSize>
 656       // from the point of parameter <param>.
 657       GCPnts_AbscissaPoint Discret( theC3d, eltSize, param );
 658       if ( !Discret.IsDone() ) break;
 659       param = Discret.Parameter();
 660       if ( param > f && param < l )
 661         theParams.push_back( param );
 662       else
 663         break;
 664       eltSize *= q;
 665     }
 666     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 667     return true;
 668   }
 669
 670   case ARITHMETIC_1D: {
 671
 672     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = theLength
 673
 674     double a1 = _value[ BEG_LENGTH_IND ];
 675     double an = _value[ END_LENGTH_IND ];
 676
 677     double  q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
 678     int     n = int( 1 + ( an - a1 ) / q );
 679
 680     double U1 = theReverse ? l : f;
 681     double Un = theReverse ? f : l;
 682     double param = U1;
 683     double eltSize = a1;
 684     if ( theReverse ) {
 685       eltSize = -eltSize;
 686       q = -q;
 687     }
 688     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 689       // computes a point on a curve <theC3d> at the distance <eltSize>
 690       // from the point of parameter <param>.
 691       GCPnts_AbscissaPoint Discret( theC3d, eltSize, param );
 692       if ( !Discret.IsDone() ) break;
 693       param = Discret.Parameter();
 694       if ( param > f && param < l )
 695         theParams.push_back( param );
 696       else
 697         break;
 698       eltSize += q;
 699     }
 700     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 701
 702     return true;
 703   }
 704
 705   case DEFLECTION: {
 706
 707     GCPnts_UniformDeflection Discret(theC3d, _value[ DEFLECTION_IND ], f, l, true);
 708     if ( !Discret.IsDone() )
 709       return false;
 710
 711     int NbPoints = Discret.NbPoints();
 712     for ( int i = 2; i < NbPoints; i++ )
 713     {
 714       double param = Discret.Parameter(i);
 715       theParams.push_back( param );
 716     }
 717     return true;
 718
 719   }
 720
 721   default:;
 722   }
 723
 724   return false;
 725 }
 726
 727 //=============================================================================
 728 /*!
 729  *
 730  */
 731 //=============================================================================
 732
 733 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
 734 {
 735   if ( _hypType == NONE )
 736     return false;
 737
 738   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 739
 740   const TopoDS_Edge & EE = TopoDS::Edge(aShape);
 741   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
 742   int shapeID = meshDS->ShapeToIndex( E );
 743
 744   double f, l;
 745   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
 746
 747   TopoDS_Vertex VFirst, VLast;
 748   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
 749
 750   ASSERT(!VFirst.IsNull());
 751   ASSERT(!VLast.IsNull());
 752   const SMDS_MeshNode * idFirst = SMESH_Algo::VertexNode( VFirst, meshDS );
 753   const SMDS_MeshNode * idLast = SMESH_Algo::VertexNode( VLast, meshDS );
 754   if (!idFirst || !idLast)
 755     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
 756
 757   if (!Curve.IsNull())
 758   {
 759     list< double > params;
 760     bool reversed = false;
 761     if ( !_mainEdge.IsNull() )
 762       reversed = aMesh.IsReversedInChain( EE, _mainEdge );
 763
 764     BRepAdaptor_Curve C3d( E );
 765     double length = EdgeLength( E );
 766     if ( ! computeInternalParameters( C3d, length, f, l, params, reversed )) {
 767       return false;
 768     }
 769     redistributeNearVertices( aMesh, C3d, length, params, VFirst, VLast );
 770
 771     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
 772     // only internal nodes receive an edge position with param on curve
 773
 774     const SMDS_MeshNode * idPrev = idFirst;
 775     double parPrev = f;
 776     double parLast = l;
 777
 778     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
 779       double param = *itU;
 780       gp_Pnt P = Curve->Value(param);
 781
 782       //Add the Node in the DataStructure
 783       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 784       meshDS->SetNodeOnEdge(node, shapeID, param);
 785
 786       if(_quadraticMesh) {
 787         // create medium node
 788         double prm = ( parPrev + param )/2;
 789         gp_Pnt PM = Curve->Value(prm);
 790         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 791         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 792         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 793         meshDS->SetMeshElementOnShape(edge, shapeID);
 794       }
 795       else {
 796         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 797         meshDS->SetMeshElementOnShape(edge, shapeID);
 798       }
 799
 800       idPrev = node;
 801       parPrev = param;
 802     }
 803     if(_quadraticMesh) {
 804       double prm = ( parPrev + parLast )/2;
 805       gp_Pnt PM = Curve->Value(prm);
 806       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 807       meshDS->SetNodeOnEdge(NM, shapeID, prm);
 808       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
 809       meshDS->SetMeshElementOnShape(edge, shapeID);
 810     }
 811     else {
 812       SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
 813       meshDS->SetMeshElementOnShape(edge, shapeID);
 814     }
 815   }
 816   else {
 817     // Edge is a degenerated Edge : We put n = 5 points on the edge.
 818     const int NbPoints = 5;
 819     BRep_Tool::Range( E, f, l ); // PAL15185
 820     double du = (l - f) / (NbPoints - 1);
 821     //MESSAGE("************* Degenerated edge! *****************");
 822
 823     gp_Pnt P = BRep_Tool::Pnt(VFirst);
 824
 825     const SMDS_MeshNode * idPrev = idFirst;
 826     for (int i = 2; i < NbPoints; i++) {
 827       double param = f + (i - 1) * du;
 828       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 829       if(_quadraticMesh) {
 830         // create medium node
 831         double prm = param - du/2.;
 832         SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
 833         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 834         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 835         meshDS->SetMeshElementOnShape(edge, shapeID);
 836       }
 837       else {
 838         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 839         meshDS->SetMeshElementOnShape(edge, shapeID);
 840       }
 841       meshDS->SetNodeOnEdge(node, shapeID, param);
 842       idPrev = node;
 843     }
 844     if(_quadraticMesh) {
 845       // create medium node
 846       double prm = l - du/2.;
 847       SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
 848       meshDS->SetNodeOnEdge(NM, shapeID, prm);
 849       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
 850       meshDS->SetMeshElementOnShape(edge, shapeID);
 851     }
 852     else {
 853       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
 854       meshDS->SetMeshElementOnShape(edge, shapeID);
 855     }
 856   }
 857   return true;
 858 }
 859
 860 //=============================================================================
 861 /*!
 862  *  See comments in SMESH_Algo.cxx
 863  */
 864 //=============================================================================
 865
 866 const list <const SMESHDS_Hypothesis *> &
 867 StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
 868                                          const TopoDS_Shape & aShape,
 869                                          const bool           ignoreAuxiliary)
 870 {
 871   _usedHypList.clear();
 872   _mainEdge.Nullify();
 873
 874   SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
 875   auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
 876   const bool ignoreAux = true;
 877   InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
 878
 879   // get non-auxiliary assigned to aShape
 880   int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
 881
 882   if (nbHyp == 0)
 883   {
 884     // Check, if propagated from some other edge
 885     if (aShape.ShapeType() == TopAbs_EDGE &&
 886         aMesh.IsPropagatedHypothesis(aShape, _mainEdge))
 887     {
 888       // Propagation of 1D hypothesis from <aMainEdge> on this edge;
 889       // get non-auxiliary assigned to _mainEdge
 890       nbHyp = aMesh.GetHypotheses( _mainEdge, compatibleFilter, _usedHypList, true );
 891     }
 892   }
 893
 894   if (nbHyp == 0) // nothing propagated nor assigned to aShape
 895   {
 896     SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
 897     nbHyp = _usedHypList.size();
 898   }
 899   else
 900   {
 901     // get auxiliary hyps from aShape
 902     aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
 903   }
 904   if ( nbHyp > 1 && ignoreAuxiliary )
 905     _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
 906
 907   return _usedHypList;
 908 }