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