src/StdMeshers/StdMeshers_Regular_1D.cxx

   1 // Copyright (C) 2007-2022  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, or (at your option) any later version.
  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 //  File   : StdMeshers_Regular_1D.cxx
  24 //           Moved here from SMESH_Regular_1D.cxx
  25 //  Author : Paul RASCLE, EDF
  26 //  Module : SMESH
  27 //
  28 #include "StdMeshers_Regular_1D.hxx"
  29
  30 #include "SMDS_MeshElement.hxx"
  31 #include "SMDS_MeshNode.hxx"
  32 #include "SMESHDS_Mesh.hxx"
  33 #include "SMESH_Comment.hxx"
  34 #include "SMESH_Gen.hxx"
  35 #include "SMESH_HypoFilter.hxx"
  36 #include "SMESH_Mesh.hxx"
  37 #include "SMESH_subMesh.hxx"
  38 #include "SMESH_subMeshEventListener.hxx"
  39 #include "StdMeshers_Adaptive1D.hxx"
  40 #include "StdMeshers_Arithmetic1D.hxx"
  41 #include "StdMeshers_AutomaticLength.hxx"
  42 #include "StdMeshers_Geometric1D.hxx"
  43 #include "StdMeshers_Deflection1D.hxx"
  44 #include "StdMeshers_Distribution.hxx"
  45 #include "StdMeshers_FixedPoints1D.hxx"
  46 #include "StdMeshers_LocalLength.hxx"
  47 #include "StdMeshers_MaxLength.hxx"
  48 #include "StdMeshers_NumberOfSegments.hxx"
  49 #include "StdMeshers_Propagation.hxx"
  50 #include "StdMeshers_SegmentLengthAroundVertex.hxx"
  51 #include "StdMeshers_StartEndLength.hxx"
  52
  53 #include <Utils_SALOME_Exception.hxx>
  54 #include <utilities.h>
  55
  56 #include <BRepAdaptor_Curve.hxx>
  57 #include <BRep_Tool.hxx>
  58 #include <GCPnts_AbscissaPoint.hxx>
  59 #include <GCPnts_UniformAbscissa.hxx>
  60 #include <GCPnts_UniformDeflection.hxx>
  61 #include <Precision.hxx>
  62 #include <TopExp.hxx>
  63 #include <TopExp_Explorer.hxx>
  64 #include <TopoDS.hxx>
  65 #include <TopoDS_Edge.hxx>
  66 #include <TopoDS_Vertex.hxx>
  67
  68 #include <string>
  69 #include <limits>
  70
  71 using namespace std;
  72 using namespace StdMeshers;
  73
  74 //=============================================================================
  75 /*!
  76  *
  77  */
  78 //=============================================================================
  79
  80 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int         hypId,
  81                                              SMESH_Gen * gen)
  82   :SMESH_1D_Algo( hypId, gen )
  83 {
  84   _name = "Regular_1D";
  85   _shapeType = (1 << TopAbs_EDGE);
  86   _fpHyp = 0;
  87
  88   _compatibleHypothesis.push_back("LocalLength");
  89   _compatibleHypothesis.push_back("MaxLength");
  90   _compatibleHypothesis.push_back("NumberOfSegments");
  91   _compatibleHypothesis.push_back("StartEndLength");
  92   _compatibleHypothesis.push_back("Deflection1D");
  93   _compatibleHypothesis.push_back("Arithmetic1D");
  94   _compatibleHypothesis.push_back("GeometricProgression");
  95   _compatibleHypothesis.push_back("FixedPoints1D");
  96   _compatibleHypothesis.push_back("AutomaticLength");
  97   _compatibleHypothesis.push_back("Adaptive1D");
  98   // auxiliary:
  99   _compatibleHypothesis.push_back("QuadraticMesh");
 100   _compatibleHypothesis.push_back("Propagation");
 101   _compatibleHypothesis.push_back("PropagOfDistribution");
 102 }
 103
 104 //=============================================================================
 105 /*!
 106  *
 107  */
 108 //=============================================================================
 109
 110 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
 111 {
 112 }
 113
 114 //=============================================================================
 115 /*!
 116  *
 117  */
 118 //=============================================================================
 119
 120 bool StdMeshers_Regular_1D::CheckHypothesis( SMESH_Mesh&         aMesh,
 121                                              const TopoDS_Shape& aShape,
 122                                              Hypothesis_Status&  aStatus )
 123 {
 124   _hypType        = NONE;
 125   _quadraticMesh  = false;
 126   _onlyUnaryInput = true;
 127
 128   // check propagation in a redefined GetUsedHypothesis()
 129   const list <const SMESHDS_Hypothesis * > & hyps =
 130     GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliaryHyps=*/false);
 131
 132   const SMESH_HypoFilter & propagFilter = StdMeshers_Propagation::GetFilter();
 133
 134   // find non-auxiliary hypothesis
 135   const SMESHDS_Hypothesis *theHyp = 0;
 136   set< string > propagTypes;
 137   list <const SMESHDS_Hypothesis * >::const_iterator h = hyps.begin();
 138   for ( ; h != hyps.end(); ++h ) {
 139     if ( static_cast<const SMESH_Hypothesis*>(*h)->IsAuxiliary() ) {
 140       if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 )
 141         _quadraticMesh = true;
 142       if ( propagFilter.IsOk( static_cast< const SMESH_Hypothesis*>( *h ), aShape ))
 143         propagTypes.insert( (*h)->GetName() );
 144     }
 145     else {
 146       if ( !theHyp )
 147         theHyp = *h; // use only the first non-auxiliary hypothesis
 148     }
 149   }
 150
 151   if ( !theHyp )
 152   {
 153     aStatus = SMESH_Hypothesis::HYP_MISSING;
 154     return false;  // can't work without a hypothesis
 155   }
 156
 157   string hypName = theHyp->GetName();
 158
 159   if ( !_mainEdge.IsNull() && _hypType == DISTRIB_PROPAGATION )
 160   {
 161     aStatus = SMESH_Hypothesis::HYP_OK;
 162   }
 163   else if ( hypName == "LocalLength" )
 164   {
 165     const StdMeshers_LocalLength * hyp =
 166       dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
 167     ASSERT(hyp);
 168     _value[ BEG_LENGTH_IND ] = hyp->GetLength();
 169     _value[ PRECISION_IND  ] = hyp->GetPrecision();
 170     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 171     _hypType = LOCAL_LENGTH;
 172     aStatus = SMESH_Hypothesis::HYP_OK;
 173   }
 174
 175   else if ( hypName == "MaxLength" )
 176   {
 177     const StdMeshers_MaxLength * hyp =
 178       dynamic_cast <const StdMeshers_MaxLength * >(theHyp);
 179     ASSERT(hyp);
 180     _value[ BEG_LENGTH_IND ] = hyp->GetLength();
 181     if ( hyp->GetUsePreestimatedLength() ) {
 182       if ( int nbSeg = aMesh.GetGen()->GetBoundaryBoxSegmentation() )
 183         _value[ BEG_LENGTH_IND ] = aMesh.GetShapeDiagonalSize() / nbSeg;
 184     }
 185     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 186     _hypType = MAX_LENGTH;
 187     aStatus = SMESH_Hypothesis::HYP_OK;
 188   }
 189
 190   else if ( hypName == "NumberOfSegments" )
 191   {
 192     const StdMeshers_NumberOfSegments * hyp =
 193       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
 194     ASSERT(hyp);
 195     _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
 196     ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
 197     _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
 198     switch (_ivalue[ DISTR_TYPE_IND ])
 199     {
 200     case StdMeshers_NumberOfSegments::DT_Scale:
 201       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
 202       _revEdgesIDs = hyp->GetReversedEdges();
 203       break;
 204     case StdMeshers_NumberOfSegments::DT_TabFunc:
 205       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
 206       _revEdgesIDs = hyp->GetReversedEdges();
 207       break;
 208     case StdMeshers_NumberOfSegments::DT_ExprFunc:
 209       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
 210       _revEdgesIDs = hyp->GetReversedEdges();
 211       break;
 212     case StdMeshers_NumberOfSegments::DT_Regular:
 213       break;
 214     default:
 215       ASSERT(0);
 216       break;
 217     }
 218     if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
 219         _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
 220         _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
 221     _hypType = NB_SEGMENTS;
 222     aStatus = SMESH_Hypothesis::HYP_OK;
 223   }
 224
 225   else if ( hypName == "Arithmetic1D" )
 226   {
 227     const StdMeshers_Arithmetic1D * hyp =
 228       dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
 229     ASSERT(hyp);
 230     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 231     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 232     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 233     _hypType = ARITHMETIC_1D;
 234
 235     _revEdgesIDs = hyp->GetReversedEdges();
 236
 237     aStatus = SMESH_Hypothesis::HYP_OK;
 238   }
 239
 240   else if ( hypName == "GeometricProgression" )
 241   {
 242     const StdMeshers_Geometric1D * hyp =
 243       dynamic_cast <const StdMeshers_Geometric1D * >(theHyp);
 244     ASSERT(hyp);
 245     _value[ BEG_LENGTH_IND ] = hyp->GetStartLength();
 246     _value[ END_LENGTH_IND ] = hyp->GetCommonRatio();
 247     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 248     _hypType = GEOMETRIC_1D;
 249
 250     _revEdgesIDs = hyp->GetReversedEdges();
 251
 252     aStatus = SMESH_Hypothesis::HYP_OK;
 253   }
 254
 255   else if ( hypName == "FixedPoints1D" ) {
 256     _fpHyp = dynamic_cast <const StdMeshers_FixedPoints1D*>(theHyp);
 257     ASSERT(_fpHyp);
 258     _hypType = FIXED_POINTS_1D;
 259
 260     _revEdgesIDs = _fpHyp->GetReversedEdges();
 261
 262     aStatus = SMESH_Hypothesis::HYP_OK;
 263   }
 264
 265   else if ( hypName == "StartEndLength" )
 266   {
 267     const StdMeshers_StartEndLength * hyp =
 268       dynamic_cast <const StdMeshers_StartEndLength * >(theHyp);
 269     ASSERT(hyp);
 270     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 271     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 272     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 273     _hypType = BEG_END_LENGTH;
 274
 275     _revEdgesIDs = hyp->GetReversedEdges();
 276
 277     aStatus = SMESH_Hypothesis::HYP_OK;
 278   }
 279
 280   else if ( hypName == "Deflection1D" )
 281   {
 282     const StdMeshers_Deflection1D * hyp =
 283       dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
 284     ASSERT(hyp);
 285     _value[ DEFLECTION_IND ] = hyp->GetDeflection();
 286     ASSERT( _value[ DEFLECTION_IND ] > 0 );
 287     _hypType = DEFLECTION;
 288     aStatus = SMESH_Hypothesis::HYP_OK;
 289   }
 290
 291   else if ( hypName == "AutomaticLength" )
 292   {
 293     StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
 294       (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
 295     ASSERT(hyp);
 296     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 297     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 298     _hypType = MAX_LENGTH;
 299     aStatus = SMESH_Hypothesis::HYP_OK;
 300   }
 301   else if ( hypName == "Adaptive1D" )
 302   {
 303     _adaptiveHyp = dynamic_cast < const StdMeshers_Adaptive1D* >(theHyp);
 304     ASSERT(_adaptiveHyp);
 305     _hypType = ADAPTIVE;
 306     _onlyUnaryInput = false;
 307     aStatus = SMESH_Hypothesis::HYP_OK;
 308   }
 309   else
 310   {
 311     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 312   }
 313
 314   if ( propagTypes.size() > 1 && aStatus == HYP_OK )
 315   {
 316     // detect concurrent Propagation hyps
 317     _usedHypList.clear();
 318     list< TopoDS_Shape > assignedTo;
 319     if ( aMesh.GetHypotheses( aShape, propagFilter, _usedHypList, true, &assignedTo ) > 1 )
 320     {
 321       // find most simple shape and a hyp on it
 322       int simpleShape = TopAbs_COMPOUND;
 323       const SMESHDS_Hypothesis* localHyp = 0;
 324       list< TopoDS_Shape >::iterator            shape = assignedTo.begin();
 325       list< const SMESHDS_Hypothesis *>::iterator hyp = _usedHypList.begin();
 326       for ( ; shape != assignedTo.end(); ++shape )
 327         if ( shape->ShapeType() > simpleShape )
 328         {
 329           simpleShape = shape->ShapeType();
 330           localHyp = (*hyp);
 331         }
 332       // check if there a different hyp on simpleShape
 333       shape = assignedTo.begin();
 334       hyp = _usedHypList.begin();
 335       for ( ; hyp != _usedHypList.end(); ++hyp, ++shape )
 336         if ( shape->ShapeType() == simpleShape &&
 337              !localHyp->IsSameName( **hyp ))
 338         {
 339           aStatus = HYP_INCOMPAT_HYPS;
 340           return error( SMESH_Comment("Hypotheses of both \"")
 341                         << StdMeshers_Propagation::GetName() << "\" and \""
 342                         << StdMeshers_PropagOfDistribution::GetName()
 343                         << "\" types can't be applied to the same edge");
 344         }
 345     }
 346   }
 347
 348   return ( aStatus == SMESH_Hypothesis::HYP_OK );
 349 }
 350
 351 static bool computeParamByFunc(Adaptor3d_Curve& C3d,
 352                                double first, double last, double length,
 353                                bool theReverse, smIdType nbSeg, Function& func,
 354                                list<double>& theParams)
 355 {
 356   // never do this way
 357   //OSD::SetSignal( true );
 358
 359   if ( nbSeg <= 0 )
 360     return false;
 361
 362   smIdType nbPnt = 1 + nbSeg;
 363   vector<double> x( nbPnt, 0. );
 364
 365
 366   const double eps = Min( 1E-4, 0.01 / double( nbSeg ));
 367
 368   if ( !buildDistribution( func, 0.0, 1.0, nbSeg, x, eps ))
 369      return false;
 370
 371   // apply parameters in range [0,1] to the space of the curve
 372   double prevU = first;
 373   double  sign = 1.;
 374   if ( theReverse )
 375   {
 376     prevU = last;
 377     sign  = -1.;
 378   }
 379
 380   for ( smIdType i = 1; i < nbSeg; i++ )
 381   {
 382     double curvLength = length * (x[i] - x[i-1]) * sign;
 383     double tol        = Min( Precision::Confusion(), curvLength / 100. );
 384     GCPnts_AbscissaPoint Discret( tol, C3d, curvLength, prevU );
 385     if ( !Discret.IsDone() )
 386       return false;
 387     double U = Discret.Parameter();
 388     if ( U > first && U < last )
 389       theParams.push_back( U );
 390     else
 391       return false;
 392     prevU = U;
 393   }
 394   if ( theReverse )
 395     theParams.reverse();
 396
 397   return true;
 398 }
 399
 400
 401 //================================================================================
 402 /*!
 403  * \brief adjust internal node parameters so that the last segment length == an
 404  *  \param a1 - the first segment length
 405  *  \param an - the last segment length
 406  *  \param U1 - the first edge parameter
 407  *  \param Un - the last edge parameter
 408  *  \param length - the edge length
 409  *  \param C3d - the edge curve
 410  *  \param theParams - internal node parameters to adjust
 411  *  \param adjustNeighbors2an - to adjust length of segments next to the last one
 412  *   and not to remove parameters
 413  */
 414 //================================================================================
 415
 416 static void compensateError(double a1, double an,
 417                             double U1, double Un,
 418                             double            length,
 419                             Adaptor3d_Curve&  C3d,
 420                             list<double> &    theParams,
 421                             bool              adjustNeighbors2an = false)
 422 {
 423   smIdType i, nPar = theParams.size();
 424   if ( a1 + an <= length && nPar > 1 )
 425   {
 426     bool reverse = ( U1 > Un );
 427     double tol   = Min( Precision::Confusion(), 0.01 * an );
 428     GCPnts_AbscissaPoint Discret( tol, C3d, reverse ? an : -an, Un );
 429     if ( !Discret.IsDone() )
 430       return;
 431     double Utgt = Discret.Parameter(); // target value of the last parameter
 432     list<double>::reverse_iterator itU = theParams.rbegin();
 433     double Ul = *itU++; // real value of the last parameter
 434     double dUn = Utgt - Ul; // parametric error of <an>
 435     double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
 436     if ( Abs(dUn) <= 1e-3 * dU )
 437       return;
 438     if ( adjustNeighbors2an || Abs(dUn) < 0.5 * dU ) { // last segment is a bit shorter than it should
 439       // move the last parameter to the edge beginning
 440     }
 441     else {  // last segment is much shorter than it should -> remove the last param and
 442       theParams.pop_back(); nPar--; // move the rest points toward the edge end
 443       dUn = Utgt - theParams.back();
 444     }
 445
 446     if ( !adjustNeighbors2an )
 447     {
 448       double q = dUn / ( Utgt - Un ); // (signed) factor of segment length change
 449       for ( itU = theParams.rbegin(), i = 1; i < nPar; i++ ) {
 450         double prevU = *itU;
 451         (*itU) += dUn;
 452         ++itU;
 453         dUn = q * (*itU - prevU) * (prevU-U1)/(Un-U1);
 454       }
 455     }
 456     else if ( nPar == 1 )
 457     {
 458       theParams.back() += dUn;
 459     }
 460     else
 461     {
 462       double q  =  dUn / double( nPar - 1 );
 463       theParams.back() += dUn;
 464       double sign = reverse ? -1 : 1;
 465       double prevU = theParams.back();
 466       itU = theParams.rbegin();
 467       for ( ++itU, i = 2; i < nPar; ++itU, i++ ) {
 468         double newU = *itU + dUn;
 469         if ( newU*sign < prevU*sign ) {
 470           prevU = *itU = newU;
 471           dUn -= q;
 472         }
 473         else { // set U between prevU and next valid param
 474           list<double>::reverse_iterator itU2 = itU;
 475           ++itU2;
 476           int nb = 2;
 477           while ( (*itU2)*sign > prevU*sign ) {
 478             ++itU2; ++nb;
 479           }
 480           dU = ( *itU2 - prevU ) / nb;
 481           while ( itU != itU2 ) {
 482             *itU += dU; ++itU;
 483           }
 484           break;
 485         }
 486       }
 487     }
 488   }
 489 }
 490
 491 //================================================================================
 492 /*!
 493  * \brief Class used to clean mesh on edges when 0D hyp modified.
 494  * Common approach doesn't work when 0D algo is missing because the 0D hyp is
 495  * considered as not participating in computation whereas it is used by 1D algo.
 496  */
 497 //================================================================================
 498
 499 // struct VertexEventListener : public SMESH_subMeshEventListener
 500 // {
 501 //   VertexEventListener():SMESH_subMeshEventListener(0) // won't be deleted by submesh
 502 //   {}
 503 //   /*!
 504 //    * \brief Clean mesh on edges
 505 //    * \param event - algo_event or compute_event itself (of SMESH_subMesh)
 506 //    * \param eventType - ALGO_EVENT or COMPUTE_EVENT (of SMESH_subMesh)
 507 //    * \param subMesh - the submesh where the event occurs
 508 //    */
 509 //   void ProcessEvent(const int event, const int eventType, SMESH_subMesh* subMesh,
 510 //                     EventListenerData*, const SMESH_Hypothesis*)
 511 //   {
 512 //     if ( eventType == SMESH_subMesh::ALGO_EVENT) // all algo events
 513 //     {
 514 //       subMesh->ComputeStateEngine( SMESH_subMesh::MODIF_ALGO_STATE );
 515 //     }
 516 //   }
 517 // }; // struct VertexEventListener
 518
 519 //=============================================================================
 520 /*!
 521  * \brief Sets event listener to vertex submeshes
 522  * \param subMesh - submesh where algo is set
 523  *
 524  * This method is called when a submesh gets HYP_OK algo_state.
 525  * After being set, event listener is notified on each event of a submesh.
 526  */
 527 //=============================================================================
 528
 529 void StdMeshers_Regular_1D::SetEventListener(SMESH_subMesh* subMesh)
 530 {
 531   StdMeshers_Propagation::SetPropagationMgr( subMesh );
 532 }
 533
 534 //=============================================================================
 535 /*!
 536  * \brief Do nothing
 537  * \param subMesh - restored submesh
 538  *
 539  * This method is called only if a submesh has HYP_OK algo_state.
 540  */
 541 //=============================================================================
 542
 543 void StdMeshers_Regular_1D::SubmeshRestored(SMESH_subMesh* /*subMesh*/)
 544 {
 545 }
 546
 547 //=============================================================================
 548 /*!
 549  * \brief Return StdMeshers_SegmentLengthAroundVertex assigned to vertex
 550  */
 551 //=============================================================================
 552
 553 const StdMeshers_SegmentLengthAroundVertex*
 554 StdMeshers_Regular_1D::getVertexHyp(SMESH_Mesh &          theMesh,
 555                                     const TopoDS_Vertex & theV)
 556 {
 557   static SMESH_HypoFilter filter( SMESH_HypoFilter::HasName("SegmentAroundVertex_0D"));
 558   if ( const SMESH_Hypothesis * h = theMesh.GetHypothesis( theV, filter, true ))
 559   {
 560     SMESH_Algo* algo = const_cast< SMESH_Algo* >( static_cast< const SMESH_Algo* > ( h ));
 561     const list <const SMESHDS_Hypothesis *> & hypList = algo->GetUsedHypothesis( theMesh, theV, 0 );
 562     if ( !hypList.empty() && string("SegmentLengthAroundVertex") == hypList.front()->GetName() )
 563       return static_cast<const StdMeshers_SegmentLengthAroundVertex*>( hypList.front() );
 564   }
 565   return 0;
 566 }
 567
 568 //================================================================================
 569 /*!
 570  * \brief Divide a curve into equal segments
 571  */
 572 //================================================================================
 573
 574 bool StdMeshers_Regular_1D::divideIntoEqualSegments( SMESH_Mesh &        theMesh,
 575                                                      Adaptor3d_Curve &   theC3d,
 576                                                      smIdType            theNbPoints,
 577                                                      double              theTol,
 578                                                      double              theLength,
 579                                                      double              theFirstU,
 580                                                      double              theLastU,
 581                                                      std::list<double> & theParameters )
 582 {
 583   bool ok = false;
 584   if ( theNbPoints < IntegerLast() )
 585   {
 586     int nbPnt = FromSmIdType<int>( theNbPoints );
 587     GCPnts_UniformAbscissa discret(theC3d, nbPnt, theFirstU, theLastU, theTol );
 588     if ( !discret.IsDone() )
 589       return error( "GCPnts_UniformAbscissa failed");
 590     if ( discret.NbPoints() < nbPnt )
 591       discret.Initialize(theC3d, nbPnt + 1, theFirstU, theLastU, theTol );
 592
 593     int nbPoints = Min( discret.NbPoints(), nbPnt );
 594     for ( int i = 2; i < nbPoints; i++ ) // skip 1st and last points
 595     {
 596       double param = discret.Parameter(i);
 597       theParameters.push_back( param );
 598     }
 599     ok = true;
 600   }
 601   else // huge nb segments
 602   {
 603     // use FIXED_POINTS_1D method
 604     StdMeshers_FixedPoints1D fixedPointsHyp( GetGen()->GetANewId(), GetGen() );
 605     _fpHyp = &fixedPointsHyp;
 606     std::vector<double>   params = { 0., 1. };
 607     std::vector<smIdType> nbSegs = { theNbPoints - 1 };
 608     fixedPointsHyp.SetPoints( params );
 609     fixedPointsHyp.SetNbSegments( nbSegs );
 610
 611     HypothesisType curType = _hypType;
 612     _hypType = FIXED_POINTS_1D;
 613
 614     ok = computeInternalParameters( theMesh, theC3d, theLength, theFirstU, theLastU,
 615                                     theParameters, /*reverse=*/false );
 616     _hypType = curType;
 617     _fpHyp = 0;
 618   }
 619   return ok;
 620 }
 621
 622 //================================================================================
 623 /*!
 624  * \brief Tune parameters to fit "SegmentLengthAroundVertex" hypothesis
 625  *  \param theC3d - wire curve
 626  *  \param theLength - curve length
 627  *  \param theParameters - internal nodes parameters to modify
 628  *  \param theVf - 1st vertex
 629  *  \param theVl - 2nd vertex
 630  */
 631 //================================================================================
 632
 633 void StdMeshers_Regular_1D::redistributeNearVertices (SMESH_Mesh &          theMesh,
 634                                                       Adaptor3d_Curve &     theC3d,
 635                                                       double                theLength,
 636                                                       std::list< double > & theParameters,
 637                                                       const TopoDS_Vertex & theVf,
 638                                                       const TopoDS_Vertex & theVl)
 639 {
 640   double f = theC3d.FirstParameter(), l = theC3d.LastParameter();
 641   size_t nPar = theParameters.size();
 642   for ( int isEnd1 = 0; isEnd1 < 2; ++isEnd1 )
 643   {
 644     const TopoDS_Vertex & V = isEnd1 ? theVf : theVl;
 645     const StdMeshers_SegmentLengthAroundVertex* hyp = getVertexHyp (theMesh, V );
 646     if ( hyp ) {
 647       double vertexLength = hyp->GetLength();
 648       if ( vertexLength > theLength / 2.0 )
 649         continue;
 650       if ( isEnd1 ) { // to have a segment of interest at end of theParameters
 651         theParameters.reverse();
 652         std::swap( f, l );
 653       }
 654       if ( _hypType == NB_SEGMENTS )
 655       {
 656         compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 657       }
 658       else if ( nPar <= 3 )
 659       {
 660         if ( !isEnd1 )
 661           vertexLength = -vertexLength;
 662         double tol = Min( Precision::Confusion(), 0.01 * vertexLength );
 663         GCPnts_AbscissaPoint Discret( tol, theC3d, vertexLength, l );
 664         if ( Discret.IsDone() ) {
 665           if ( nPar == 0 )
 666             theParameters.push_back( Discret.Parameter());
 667           else {
 668             double L = GCPnts_AbscissaPoint::Length( theC3d, theParameters.back(), l);
 669             if ( vertexLength < L / 2.0 )
 670               theParameters.push_back( Discret.Parameter());
 671             else
 672               compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 673           }
 674         }
 675       }
 676       else
 677       {
 678         // recompute params between the last segment and a middle one.
 679         // find size of a middle segment
 680         smIdType nHalf = ( nPar-1 ) / 2;
 681         list< double >::reverse_iterator itU = theParameters.rbegin();
 682         std::advance( itU, nHalf );
 683         double Um = *itU++;
 684         double Lm = GCPnts_AbscissaPoint::Length( theC3d, Um, *itU);
 685         double L = GCPnts_AbscissaPoint::Length( theC3d, *itU, l);
 686         static StdMeshers_Regular_1D* auxAlgo = 0;
 687         if ( !auxAlgo ) {
 688           auxAlgo = new StdMeshers_Regular_1D( _gen->GetANewId(), _gen );
 689           auxAlgo->_hypType = BEG_END_LENGTH;
 690         }
 691         auxAlgo->_value[ BEG_LENGTH_IND ] = Lm;
 692         auxAlgo->_value[ END_LENGTH_IND ] = vertexLength;
 693         double from = *itU, to = l;
 694         if ( isEnd1 ) {
 695           std::swap( from, to );
 696           std::swap( auxAlgo->_value[ BEG_LENGTH_IND ], auxAlgo->_value[ END_LENGTH_IND ]);
 697         }
 698         list<double> params;
 699         if ( auxAlgo->computeInternalParameters( theMesh, theC3d, L, from, to, params, false ))
 700         {
 701           if ( isEnd1 ) params.reverse();
 702           while ( 1 + nHalf-- )
 703             theParameters.pop_back();
 704           theParameters.splice( theParameters.end(), params );
 705         }
 706         else
 707         {
 708           compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 709         }
 710       }
 711       if ( isEnd1 )
 712         theParameters.reverse();
 713     }
 714   }
 715 }
 716
 717 //=============================================================================
 718 /*!
 719  *
 720  */
 721 //=============================================================================
 722 bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 723                                                       Adaptor3d_Curve& theC3d,
 724                                                       double           theLength,
 725                                                       double           theFirstU,
 726                                                       double           theLastU,
 727                                                       list<double> &   theParams,
 728                                                       const bool       theReverse,
 729                                                       bool             theConsiderPropagation)
 730 {
 731   theParams.clear();
 732
 733   double f = theFirstU, l = theLastU;
 734
 735   // Propagation Of Distribution
 736   //
 737   if ( !_mainEdge.IsNull() && _hypType == DISTRIB_PROPAGATION )
 738   {
 739     TopoDS_Edge mainEdge = TopoDS::Edge( _mainEdge ); // should not be a reference!
 740     _gen->Compute( theMesh, mainEdge, SMESH_Gen::SHAPE_ONLY_UPWARD );
 741
 742     SMESHDS_SubMesh* smDS = theMesh.GetMeshDS()->MeshElements( mainEdge );
 743     if ( !smDS )
 744       return error("No mesh on the source edge of Propagation Of Distribution");
 745     if ( smDS->NbNodes() < 1 )
 746       return true; // 1 segment
 747
 748     map< double, const SMDS_MeshNode* > mainEdgeParamsOfNodes;
 749     if ( ! SMESH_Algo::GetSortedNodesOnEdge( theMesh.GetMeshDS(), mainEdge, _quadraticMesh,
 750                                              mainEdgeParamsOfNodes, SMDSAbs_Edge ))
 751       return error("Bad node parameters on the source edge of Propagation Of Distribution");
 752     vector< double > segLen( mainEdgeParamsOfNodes.size() - 1 );
 753     double totalLen = 0;
 754     BRepAdaptor_Curve mainEdgeCurve( mainEdge );
 755     map< double, const SMDS_MeshNode* >::iterator
 756       u_n2 = mainEdgeParamsOfNodes.begin(), u_n1 = u_n2++;
 757     for ( size_t i = 1; i < mainEdgeParamsOfNodes.size(); ++i, ++u_n1, ++u_n2 )
 758     {
 759       segLen[ i-1 ] = GCPnts_AbscissaPoint::Length( mainEdgeCurve,
 760                                                     u_n1->first,
 761                                                     u_n2->first);
 762       totalLen += segLen[ i-1 ];
 763     }
 764     for ( size_t i = 0; i < segLen.size(); ++i )
 765       segLen[ i ] *= theLength / totalLen;
 766
 767     size_t  iSeg = theReverse ? segLen.size()-1 : 0;
 768     size_t  dSeg = theReverse ? -1 : +1;
 769     double param = theFirstU;
 770     size_t nbParams = 0;
 771     for ( size_t i = 1; i < segLen.size(); ++i, iSeg += dSeg )
 772     {
 773       double tol = Min( Precision::Confusion(), 0.01 * segLen[ iSeg ]);
 774       GCPnts_AbscissaPoint Discret( tol, theC3d, segLen[ iSeg ], param );
 775       if ( !Discret.IsDone() ) break;
 776       param = Discret.Parameter();
 777       theParams.push_back( param );
 778       ++nbParams;
 779     }
 780     if ( nbParams != segLen.size()-1 )
 781       return error( SMESH_Comment("Can't divide into ") << segLen.size() << " segments");
 782
 783     compensateError( segLen[ theReverse ? segLen.size()-1 : 0 ],
 784                      segLen[ theReverse ? 0 : segLen.size()-1 ],
 785                      f, l, theLength, theC3d, theParams, true );
 786     return true;
 787   }
 788
 789
 790   switch( _hypType )
 791   {
 792   case LOCAL_LENGTH:
 793   case MAX_LENGTH:
 794   case NB_SEGMENTS:
 795   {
 796     double eltSize = 1;
 797     smIdType nbSegments;
 798     if ( _hypType == MAX_LENGTH )
 799     {
 800       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 801       if (nbseg <= 0)
 802         nbseg = 1; // degenerated edge
 803       eltSize = theLength / nbseg * ( 1. - 1e-9 );
 804       nbSegments = ToSmIdType( nbseg );
 805     }
 806     else if ( _hypType == LOCAL_LENGTH )
 807     {
 808       // Local Length hypothesis
 809       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 810
 811       // NPAL17873:
 812       bool isFound = false;
 813       if (theConsiderPropagation && !_mainEdge.IsNull()) // propagated from some other edge
 814       {
 815         // Advanced processing to assure equal number of segments in case of Propagation
 816         SMESH_subMesh* sm = theMesh.GetSubMeshContaining(_mainEdge);
 817         if (sm) {
 818           bool computed = sm->IsMeshComputed();
 819           if (!computed) {
 820             if (sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE) {
 821               _gen->Compute( theMesh, _mainEdge, /*anUpward=*/true);
 822               computed = sm->IsMeshComputed();
 823             }
 824           }
 825           if (computed) {
 826             SMESHDS_SubMesh* smds = sm->GetSubMeshDS();
 827             smIdType  nb_segments = smds->NbElements();
 828             if (nbseg - 1 <= nb_segments && nb_segments <= nbseg + 1) {
 829               isFound = true;
 830               nbseg = FromSmIdType<double>( nb_segments );
 831             }
 832           }
 833         }
 834       }
 835       if (!isFound) // not found by meshed edge in the propagation chain, use precision
 836       {
 837         double aPrecision = _value[ PRECISION_IND ];
 838         double nbseg_prec = ceil((theLength / _value[ BEG_LENGTH_IND ]) - aPrecision);
 839         if (nbseg_prec == (nbseg - 1)) nbseg--;
 840       }
 841
 842       if (nbseg <= 0)
 843         nbseg = 1;                        // degenerated edge
 844       eltSize = theLength / nbseg;
 845       nbSegments = ToSmIdType( nbseg );
 846     }
 847     else
 848     {
 849       // Number Of Segments hypothesis
 850       nbSegments = _ivalue[ NB_SEGMENTS_IND ];
 851       if ( nbSegments < 1 )  return false;
 852       if ( nbSegments == 1 ) return true;
 853
 854       switch (_ivalue[ DISTR_TYPE_IND ])
 855       {
 856       case StdMeshers_NumberOfSegments::DT_Scale:
 857         {
 858           double scale = _value[ SCALE_FACTOR_IND ];
 859
 860           if (fabs(scale - 1.0) < Precision::Confusion()) {
 861             // special case to avoid division by zero
 862             for ( smIdType i = 1; i < nbSegments; i++) {
 863               double param = f + (l - f) * double( i ) / double( nbSegments );
 864               theParams.push_back( param );
 865             }
 866           }
 867           else { // general case of scale distribution
 868             if ( theReverse )
 869               scale = 1.0 / scale;
 870
 871             double  alpha = pow(scale, 1.0 / double( nbSegments - 1 ));
 872             double factor = (l - f) / (1.0 - pow(alpha, nbSegments));
 873
 874             for ( smIdType i = 1; i < nbSegments; i++) {
 875               double param = f + factor * (1.0 - pow(alpha, i));
 876               theParams.push_back( param );
 877             }
 878           }
 879           const double lenFactor = theLength/(l-f);
 880           const double minSegLen = Min( theParams.front() - f, l - theParams.back() );
 881           const double       tol = Min( Precision::Confusion(), 0.01 * minSegLen );
 882           list<double>::iterator u = theParams.begin(), uEnd = theParams.end();
 883           for ( ; u != uEnd; ++u )
 884           {
 885             GCPnts_AbscissaPoint Discret( tol, theC3d, ((*u)-f) * lenFactor, f );
 886             if ( Discret.IsDone() )
 887               *u = Discret.Parameter();
 888           }
 889           return true;
 890         }
 891         break;
 892       case StdMeshers_NumberOfSegments::DT_TabFunc:
 893         {
 894           FunctionTable func(_vvalue[ TAB_FUNC_IND ], FromSmIdType<int>( _ivalue[ CONV_MODE_IND ]));
 895           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 896                                     _ivalue[ NB_SEGMENTS_IND ], func,
 897                                     theParams);
 898         }
 899         break;
 900       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 901         {
 902           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(),
 903                             FromSmIdType<int>( _ivalue[ CONV_MODE_IND ]));
 904           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 905                                     _ivalue[ NB_SEGMENTS_IND ], func,
 906                                     theParams);
 907         }
 908         break;
 909       case StdMeshers_NumberOfSegments::DT_Regular:
 910         eltSize = theLength / double( nbSegments );
 911         break;
 912       default:
 913         return false;
 914       }
 915     }
 916
 917     double tol = Min( Precision::Confusion(), 0.01 * eltSize );
 918     divideIntoEqualSegments( theMesh, theC3d, nbSegments + 1, tol,
 919                              theLength, theFirstU, theLastU, theParams );
 920
 921     compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams, true ); // for PAL9899
 922     return true;
 923   }
 924
 925
 926   case BEG_END_LENGTH: {
 927
 928     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = theLength
 929
 930     double a1 = _value[ BEG_LENGTH_IND ];
 931     double an = _value[ END_LENGTH_IND ];
 932     double q  = ( theLength - a1 ) / ( theLength - an );
 933     if ( q < theLength/1e6 || 1.01*theLength < a1 + an)
 934       return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
 935                      "for an edge of length "<<theLength);
 936
 937     double      U1 = theReverse ? l : f;
 938     double      Un = theReverse ? f : l;
 939     double   param = U1;
 940     double eltSize = theReverse ? -a1 : a1;
 941     double     tol = Min( Precision::Confusion(), 0.01 * Min( a1, an ));
 942     while ( 1 ) {
 943       // computes a point on a curve <theC3d> at the distance <eltSize>
 944       // from the point of parameter <param>.
 945       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
 946       if ( !Discret.IsDone() ) break;
 947       param = Discret.Parameter();
 948       if ( f < param && param < l )
 949         theParams.push_back( param );
 950       else
 951         break;
 952       eltSize *= q;
 953     }
 954     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 955     if (theReverse) theParams.reverse(); // NPAL18025
 956     return true;
 957   }
 958
 959   case ARITHMETIC_1D:
 960   {
 961     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = theLength
 962
 963     double a1 = _value[ BEG_LENGTH_IND ];
 964     double an = _value[ END_LENGTH_IND ];
 965     if ( 1.01*theLength < a1 + an )
 966       return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
 967                      "for an edge of length "<<theLength);
 968
 969     double q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
 970     int    n = int(fabs(q) > numeric_limits<double>::min() ? ( 1+( an-a1 )/q ) : ( 1+theLength/a1 ));
 971
 972     double      U1 = theReverse ? l : f;
 973     double      Un = theReverse ? f : l;
 974     double   param = U1;
 975     double eltSize = a1;
 976     double     tol = Min( Precision::Confusion(), 0.01 * Min( a1, an ));
 977     if ( theReverse ) {
 978       eltSize = -eltSize;
 979       q = -q;
 980     }
 981     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 982       // computes a point on a curve <theC3d> at the distance <eltSize>
 983       // from the point of parameter <param>.
 984       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
 985       if ( !Discret.IsDone() ) break;
 986       param = Discret.Parameter();
 987       if ( param > f && param < l )
 988         theParams.push_back( param );
 989       else
 990         break;
 991       eltSize += q;
 992     }
 993     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 994     if ( theReverse ) theParams.reverse(); // NPAL18025
 995
 996     return true;
 997   }
 998
 999   case GEOMETRIC_1D:
1000   {
1001     double a1 = _value[ BEG_LENGTH_IND ], an = 0;
1002     double q  = _value[ END_LENGTH_IND ];
1003
1004     double U1 = theReverse ? l : f;
1005     double Un = theReverse ? f : l;
1006     double param = U1;
1007     double eltSize = a1;
1008     if ( theReverse )
1009       eltSize = -eltSize;
1010
1011     int nbParams = 0;
1012     while ( true ) {
1013       // computes a point on a curve <theC3d> at the distance <eltSize>
1014       // from the point of parameter <param>.
1015       double tol = Min( Precision::Confusion(), 0.01 * eltSize );
1016       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
1017       if ( !Discret.IsDone() ) break;
1018       param = Discret.Parameter();
1019       if ( f < param && param < l )
1020         theParams.push_back( param );
1021       else
1022         break;
1023       an = eltSize;
1024       eltSize *= q;
1025       ++nbParams;
1026       if ( q < 1. && eltSize < 1e-100 )
1027         return error("Too small common ratio causes too many segments");
1028     }
1029     if ( nbParams > 1 )
1030     {
1031       if ( Abs( param - Un ) < 0.2 * Abs( param - theParams.back() ))
1032       {
1033         compensateError( a1, Abs(eltSize), U1, Un, theLength, theC3d, theParams );
1034       }
1035       else if ( Abs( Un - theParams.back() ) <
1036                 0.2 * Abs( theParams.back() - *(++theParams.rbegin())))
1037       {
1038         theParams.pop_back();
1039         compensateError( a1, Abs(an), U1, Un, theLength, theC3d, theParams );
1040       }
1041     }
1042     if (theReverse) theParams.reverse(); // NPAL18025
1043
1044     return true;
1045   }
1046
1047   case FIXED_POINTS_1D:
1048   {
1049     const std::vector<double>& aPnts = _fpHyp->GetPoints();
1050     std::vector<smIdType>     nbsegs = _fpHyp->GetNbSegments();
1051
1052     // sort normalized params, taking into account theReverse
1053     TColStd_SequenceOfReal Params;
1054     double tol = 1e-7;
1055     for ( size_t i = 0; i < aPnts.size(); i++ )
1056     {
1057       if( aPnts[i] < tol || aPnts[i] > 1 - tol )
1058         continue;
1059       double u = theReverse ? ( 1 - aPnts[i] ) : aPnts[i];
1060       int    j = 1;
1061       bool IsExist = false;
1062       for ( ; j <= Params.Length() &&  !IsExist; j++ )
1063       {
1064         IsExist = ( Abs( u - Params.Value(j) ) < tol );
1065         if ( u < Params.Value(j) ) break;
1066       }
1067       if ( !IsExist ) Params.InsertBefore( j, u );
1068     }
1069     Params.InsertBefore( 1, 0.0 );
1070     Params.Append( 1.0 );
1071
1072     if ( theReverse )
1073     {
1074       if ((int) nbsegs.size() > Params.Length() - 1 )
1075         nbsegs.resize( Params.Length() - 1 );
1076       std::reverse( nbsegs.begin(), nbsegs.end() );
1077     }
1078     if ( nbsegs.empty() )
1079     {
1080       nbsegs.push_back( 1 );
1081     }
1082     if ((int) nbsegs.size() < Params.Length() - 1 )
1083       nbsegs.resize( Params.Length() - 1, nbsegs[0] );
1084
1085     // care of huge nbsegs - additionally divide diapasons
1086     for ( int i = 2; i <= Params.Length(); i++ )
1087     {
1088       smIdType nbTot = nbsegs[ i-2 ];
1089       if ( nbTot <= IntegerLast() )
1090         continue;
1091       smIdType   nbDiapason = nbTot / IntegerLast() + 1;
1092       smIdType nbSegPerDiap = nbTot / nbDiapason;
1093       double           par0 = Params( i - 1 ), par1 = Params( i );
1094       for ( smIdType iDiap = 0; iDiap < nbDiapason - 1; ++iDiap )
1095       {
1096         double    r = double( nbSegPerDiap * ( iDiap + 1 )) / double( nbTot );
1097         double parI = par0 + ( par1 - par0 ) * r;
1098         Params.InsertBefore( i, parI );
1099         auto it = nbsegs.begin();
1100         smIdType incr_it = i - 2 + iDiap;
1101         nbsegs.insert( it + incr_it, nbSegPerDiap );
1102       }
1103       nbsegs[ i-2 + nbDiapason - 1 ] = nbSegPerDiap + nbTot % nbDiapason;
1104     }
1105
1106     // transform normalized Params into real ones
1107     std::vector< double > uVec( Params.Length() );
1108     uVec[ 0 ] = theFirstU;
1109     double abscissa;
1110     for ( int i = 2; i < Params.Length(); i++ )
1111     {
1112       abscissa = Params( i ) * theLength;
1113       tol      = Min( Precision::Confusion(), 0.01 * abscissa );
1114       GCPnts_AbscissaPoint APnt( tol, theC3d, abscissa, theFirstU );
1115       if ( !APnt.IsDone() )
1116         return error( "GCPnts_AbscissaPoint failed");
1117       uVec[ i-1 ] = APnt.Parameter();
1118     }
1119     uVec.back() = theLastU;
1120
1121     // divide segments
1122     double eltSize, segmentSize, par1, par2;
1123     for ( int i = 0; i < (int)uVec.size()-1; i++ )
1124     {
1125       par1 = uVec[ i   ];
1126       par2 = uVec[ i+1 ];
1127       smIdType nbseg = ( i < (int) nbsegs.size() ) ? nbsegs[i] : nbsegs[0];
1128       if ( nbseg > 1 )
1129       {
1130         segmentSize = ( Params( i+2 ) - Params( i+1 )) * theLength;
1131         eltSize     = segmentSize / double( nbseg );
1132         tol         = Min( Precision::Confusion(), 0.01 * eltSize );
1133         if ( !divideIntoEqualSegments( theMesh, theC3d, nbseg + 1, tol,
1134                                        segmentSize, par1, par2, theParams ))
1135           return false;
1136       }
1137       theParams.push_back( par2 );
1138     }
1139     theParams.pop_back();
1140
1141     return true;
1142   }
1143
1144   case DEFLECTION:
1145   {
1146     GCPnts_UniformDeflection Discret( theC3d, _value[ DEFLECTION_IND ], f, l, true );
1147     if ( !Discret.IsDone() )
1148       return false;
1149
1150     int NbPoints = Discret.NbPoints();
1151     for ( int i = 2; i < NbPoints; i++ )
1152     {
1153       double param = Discret.Parameter(i);
1154       theParams.push_back( param );
1155     }
1156     return true;
1157   }
1158
1159   default:;
1160   }
1161
1162   return false;
1163 }
1164
1165 //=============================================================================
1166 /*!
1167  *
1168  */
1169 //=============================================================================
1170
1171 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & theShape)
1172 {
1173   if ( _hypType == NONE )
1174     return false;
1175
1176   if ( _hypType == ADAPTIVE )
1177   {
1178     _adaptiveHyp->GetAlgo()->InitComputeError();
1179     _adaptiveHyp->GetAlgo()->Compute( theMesh, theShape );
1180     return error( _adaptiveHyp->GetAlgo()->GetComputeError() );
1181   }
1182
1183   SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
1184
1185   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
1186   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
1187   int shapeID = meshDS->ShapeToIndex( E );
1188
1189   double f, l;
1190   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
1191
1192   TopoDS_Vertex VFirst, VLast;
1193   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
1194
1195   ASSERT(!VFirst.IsNull());
1196   ASSERT(!VLast.IsNull());
1197   const SMDS_MeshNode * nFirst = SMESH_Algo::VertexNode( VFirst, meshDS );
1198   const SMDS_MeshNode *  nLast = SMESH_Algo::VertexNode( VLast,  meshDS );
1199   if ( !nFirst || !nLast )
1200     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
1201
1202   // remove elements created by e.g. pattern mapping (PAL21999)
1203   // CLEAN event is incorrectly ptopagated seemingly due to Propagation hyp
1204   // so TEMPORARY solution is to clean the submesh manually
1205   if (SMESHDS_SubMesh * subMeshDS = meshDS->MeshElements(theShape))
1206   {
1207     SMDS_ElemIteratorPtr ite = subMeshDS->GetElements();
1208     while (ite->more())
1209       meshDS->RemoveFreeElement(ite->next(), subMeshDS);
1210     SMDS_NodeIteratorPtr itn = subMeshDS->GetNodes();
1211     while (itn->more()) {
1212       const SMDS_MeshNode * node = itn->next();
1213       if ( node->NbInverseElements() == 0 )
1214         meshDS->RemoveFreeNode(node, subMeshDS);
1215       else
1216         meshDS->RemoveNode(node);
1217     }
1218   }
1219
1220   double length = EdgeLength( E );
1221   if ( !Curve.IsNull() && length > 0 )
1222   {
1223     list< double > params;
1224     bool reversed = false;
1225     if ( theMesh.GetShapeToMesh().ShapeType() >= TopAbs_WIRE && _revEdgesIDs.empty() ) {
1226       // if the shape to mesh is WIRE or EDGE
1227       reversed = ( EE.Orientation() == TopAbs_REVERSED );
1228     }
1229     if ( !_mainEdge.IsNull() ) {
1230       // take into account reversing the edge the hypothesis is propagated from
1231       // (_mainEdge.Orientation() marks mutual orientation of EDGEs in propagation chain)
1232       reversed = ( _mainEdge.Orientation() == TopAbs_REVERSED );
1233       if ( _hypType != DISTRIB_PROPAGATION ) {
1234         int mainID = meshDS->ShapeToIndex(_mainEdge);
1235         if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), mainID) != _revEdgesIDs.end())
1236           reversed = !reversed;
1237       }
1238     }
1239     // take into account this edge reversing
1240     if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), shapeID) != _revEdgesIDs.end())
1241       reversed = !reversed;
1242
1243     BRepAdaptor_Curve C3d( E );
1244     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, reversed, true )) {
1245       return false;
1246     }
1247     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
1248
1249     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
1250     // only internal nodes receive an edge position with param on curve
1251
1252     const SMDS_MeshNode * nPrev = nFirst;
1253     double parPrev = f;
1254     double parLast = l;
1255
1256     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
1257       double param = *itU;
1258       gp_Pnt P = Curve->Value(param);
1259
1260       //Add the Node in the DataStructure
1261       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
1262       meshDS->SetNodeOnEdge(node, shapeID, param);
1263
1264       if(_quadraticMesh) {
1265         // create medium node
1266         double prm = ( parPrev + param )/2;
1267         gp_Pnt  PM = Curve->Value(prm);
1268         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
1269         meshDS->SetNodeOnEdge(NM, shapeID, prm);
1270         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node, NM);
1271         meshDS->SetMeshElementOnShape(edge, shapeID);
1272       }
1273       else {
1274         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node);
1275         meshDS->SetMeshElementOnShape(edge, shapeID);
1276       }
1277
1278       nPrev   = node;
1279       parPrev = param;
1280     }
1281     if(_quadraticMesh) {
1282       double prm = ( parPrev + parLast )/2;
1283       gp_Pnt PM = Curve->Value(prm);
1284       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
1285       meshDS->SetNodeOnEdge(NM, shapeID, prm);
1286       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast, NM);
1287       meshDS->SetMeshElementOnShape(edge, shapeID);
1288     }
1289     else {
1290       SMDS_MeshEdge* edge = meshDS->AddEdge(nPrev, nLast);
1291       meshDS->SetMeshElementOnShape(edge, shapeID);
1292     }
1293   }
1294   else
1295   {
1296     // Edge is a degenerated Edge : We put n = 5 points on the edge.
1297     const int NbPoints = 5;
1298     BRep_Tool::Range( E, f, l ); // PAL15185
1299     double du = (l - f) / (NbPoints - 1);
1300
1301     gp_Pnt P = BRep_Tool::Pnt(VFirst);
1302
1303     const SMDS_MeshNode * nPrev = nFirst;
1304     for (int i = 2; i < NbPoints; i++) {
1305       double param = f + (i - 1) * du;
1306       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
1307       if(_quadraticMesh) {
1308         // create medium node
1309         double prm = param - du/2.;
1310         SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
1311         meshDS->SetNodeOnEdge(NM, shapeID, prm);
1312         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node, NM);
1313         meshDS->SetMeshElementOnShape(edge, shapeID);
1314       }
1315       else {
1316         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node);
1317         meshDS->SetMeshElementOnShape(edge, shapeID);
1318       }
1319       meshDS->SetNodeOnEdge(node, shapeID, param);
1320       nPrev = node;
1321     }
1322     if(_quadraticMesh) {
1323       // create medium node
1324       double prm = l - du/2.;
1325       SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
1326       meshDS->SetNodeOnEdge(NM, shapeID, prm);
1327       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast, NM);
1328       meshDS->SetMeshElementOnShape(edge, shapeID);
1329     }
1330     else {
1331       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast);
1332       meshDS->SetMeshElementOnShape(edge, shapeID);
1333     }
1334   }
1335   return true;
1336 }
1337
1338
1339 //=============================================================================
1340 /*!
1341  *
1342  */
1343 //=============================================================================
1344
1345 bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh &         theMesh,
1346                                      const TopoDS_Shape & theShape,
1347                                      MapShapeNbElems&     theResMap)
1348 {
1349   if ( _hypType == NONE )
1350     return false;
1351
1352   if ( _hypType == ADAPTIVE )
1353   {
1354     _adaptiveHyp->GetAlgo()->InitComputeError();
1355     _adaptiveHyp->GetAlgo()->Evaluate( theMesh, theShape, theResMap );
1356     return error( _adaptiveHyp->GetAlgo()->GetComputeError() );
1357   }
1358
1359   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
1360   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
1361
1362   double f, l;
1363   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
1364
1365   TopoDS_Vertex VFirst, VLast;
1366   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
1367
1368   ASSERT(!VFirst.IsNull());
1369   ASSERT(!VLast.IsNull());
1370
1371   std::vector<smIdType> aVec(SMDSEntity_Last,0);
1372
1373   double length = EdgeLength( E );
1374   if ( !Curve.IsNull() && length > 0 )
1375   {
1376     list< double > params;
1377     BRepAdaptor_Curve C3d( E );
1378     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, false, true )) {
1379       SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1380       theResMap.insert(std::make_pair(sm,aVec));
1381       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1382       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1383       return false;
1384     }
1385     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
1386
1387     if(_quadraticMesh) {
1388       aVec[SMDSEntity_Node     ] = 2*params.size() + 1;
1389       aVec[SMDSEntity_Quad_Edge] = params.size() + 1;
1390     }
1391     else {
1392       aVec[SMDSEntity_Node] = params.size();
1393       aVec[SMDSEntity_Edge] = params.size() + 1;
1394     }
1395
1396   }
1397   else {
1398     // Edge is a degenerated Edge : We put n = 5 points on the edge.
1399     if ( _quadraticMesh ) {
1400       aVec[SMDSEntity_Node     ] = 11;
1401       aVec[SMDSEntity_Quad_Edge] = 6;
1402     }
1403     else {
1404       aVec[SMDSEntity_Node] = 5;
1405       aVec[SMDSEntity_Edge] = 6;
1406     }
1407   }
1408
1409   SMESH_subMesh * sm = theMesh.GetSubMesh( theShape );
1410   theResMap.insert( std::make_pair( sm, aVec ));
1411
1412   return true;
1413 }
1414
1415
1416 //=============================================================================
1417 /*!
1418  *  See comments in SMESH_Algo.cxx
1419  */
1420 //=============================================================================
1421
1422 const list <const SMESHDS_Hypothesis *> &
1423 StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
1424                                          const TopoDS_Shape & aShape,
1425                                          const bool           ignoreAuxiliary)
1426 {
1427   _usedHypList.clear();
1428   _mainEdge.Nullify();
1429
1430   SMESH_HypoFilter auxiliaryFilter( SMESH_HypoFilter::IsAuxiliary() );
1431   const SMESH_HypoFilter* compatibleFilter = GetCompatibleHypoFilter(/*ignoreAux=*/true );
1432
1433   // get non-auxiliary assigned directly to aShape
1434   int nbHyp = aMesh.GetHypotheses( aShape, *compatibleFilter, _usedHypList, false );
1435
1436   if (nbHyp == 0 && aShape.ShapeType() == TopAbs_EDGE)
1437   {
1438     // Check, if propagated from some other edge
1439     bool isPropagOfDistribution = false;
1440     _mainEdge = StdMeshers_Propagation::GetPropagationSource( aMesh, aShape,
1441                                                               isPropagOfDistribution );
1442     if ( !_mainEdge.IsNull() )
1443     {
1444       if ( isPropagOfDistribution )
1445         _hypType = DISTRIB_PROPAGATION;
1446       // Propagation of 1D hypothesis from <aMainEdge> on this edge;
1447       // get non-auxiliary assigned to _mainEdge
1448       nbHyp = aMesh.GetHypotheses( _mainEdge, *compatibleFilter, _usedHypList, true );
1449     }
1450   }
1451
1452   if (nbHyp == 0) // nothing propagated nor assigned to aShape
1453   {
1454     SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
1455     nbHyp = (int)_usedHypList.size();
1456   }
1457   else
1458   {
1459     // get auxiliary hyps from aShape
1460     aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
1461   }
1462   if ( nbHyp > 1 && ignoreAuxiliary )
1463     _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
1464
1465   return _usedHypList;
1466 }
1467
1468 //================================================================================
1469 /*!
1470  * \brief Pass CancelCompute() to a child algorithm
1471  */
1472 //================================================================================
1473
1474 void StdMeshers_Regular_1D::CancelCompute()
1475 {
1476   SMESH_Algo::CancelCompute();
1477   if ( _hypType == ADAPTIVE )
1478     _adaptiveHyp->GetAlgo()->CancelCompute();
1479 }