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