src/StdMeshers/StdMeshers_Regular_1D.cxx

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