src/StdMeshers/StdMeshers_Regular_1D.cxx

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