src/StdMeshers/StdMeshers_Regular_1D.cxx

   1 // Copyright (C) 2007-2021  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, int 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   int nbPnt = 1 + nbSeg;
 363   vector<double> x( nbPnt, 0. );
 364
 365
 366   const double eps = Min( 1E-4, 1./nbSeg/100. );
 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 ( int 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   int 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 / ( 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 Tune parameters to fit "SegmentLengthAroundVertex" hypothesis
 571   * \param theC3d - wire curve
 572   * \param theLength - curve length
 573   * \param theParameters - internal nodes parameters to modify
 574   * \param theVf - 1st vertex
 575   * \param theVl - 2nd vertex
 576  */
 577 //================================================================================
 578
 579 void StdMeshers_Regular_1D::redistributeNearVertices (SMESH_Mesh &          theMesh,
 580                                                       Adaptor3d_Curve &     theC3d,
 581                                                       double                theLength,
 582                                                       std::list< double > & theParameters,
 583                                                       const TopoDS_Vertex & theVf,
 584                                                       const TopoDS_Vertex & theVl)
 585 {
 586   double f = theC3d.FirstParameter(), l = theC3d.LastParameter();
 587   int nPar = theParameters.size();
 588   for ( int isEnd1 = 0; isEnd1 < 2; ++isEnd1 )
 589   {
 590     const TopoDS_Vertex & V = isEnd1 ? theVf : theVl;
 591     const StdMeshers_SegmentLengthAroundVertex* hyp = getVertexHyp (theMesh, V );
 592     if ( hyp ) {
 593       double vertexLength = hyp->GetLength();
 594       if ( vertexLength > theLength / 2.0 )
 595         continue;
 596       if ( isEnd1 ) { // to have a segment of interest at end of theParameters
 597         theParameters.reverse();
 598         std::swap( f, l );
 599       }
 600       if ( _hypType == NB_SEGMENTS )
 601       {
 602         compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 603       }
 604       else if ( nPar <= 3 )
 605       {
 606         if ( !isEnd1 )
 607           vertexLength = -vertexLength;
 608         double tol = Min( Precision::Confusion(), 0.01 * vertexLength );
 609         GCPnts_AbscissaPoint Discret( tol, theC3d, vertexLength, l );
 610         if ( Discret.IsDone() ) {
 611           if ( nPar == 0 )
 612             theParameters.push_back( Discret.Parameter());
 613           else {
 614             double L = GCPnts_AbscissaPoint::Length( theC3d, theParameters.back(), l);
 615             if ( vertexLength < L / 2.0 )
 616               theParameters.push_back( Discret.Parameter());
 617             else
 618               compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 619           }
 620         }
 621       }
 622       else
 623       {
 624         // recompute params between the last segment and a middle one.
 625         // find size of a middle segment
 626         int nHalf = ( nPar-1 ) / 2;
 627         list< double >::reverse_iterator itU = theParameters.rbegin();
 628         std::advance( itU, nHalf );
 629         double Um = *itU++;
 630         double Lm = GCPnts_AbscissaPoint::Length( theC3d, Um, *itU);
 631         double L = GCPnts_AbscissaPoint::Length( theC3d, *itU, l);
 632         static StdMeshers_Regular_1D* auxAlgo = 0;
 633         if ( !auxAlgo ) {
 634           auxAlgo = new StdMeshers_Regular_1D( _gen->GetANewId(), _gen );
 635           auxAlgo->_hypType = BEG_END_LENGTH;
 636         }
 637         auxAlgo->_value[ BEG_LENGTH_IND ] = Lm;
 638         auxAlgo->_value[ END_LENGTH_IND ] = vertexLength;
 639         double from = *itU, to = l;
 640         if ( isEnd1 ) {
 641           std::swap( from, to );
 642           std::swap( auxAlgo->_value[ BEG_LENGTH_IND ], auxAlgo->_value[ END_LENGTH_IND ]);
 643         }
 644         list<double> params;
 645         if ( auxAlgo->computeInternalParameters( theMesh, theC3d, L, from, to, params, false ))
 646         {
 647           if ( isEnd1 ) params.reverse();
 648           while ( 1 + nHalf-- )
 649             theParameters.pop_back();
 650           theParameters.splice( theParameters.end(), params );
 651         }
 652         else
 653         {
 654           compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 655         }
 656       }
 657       if ( isEnd1 )
 658         theParameters.reverse();
 659     }
 660   }
 661 }
 662
 663 //=============================================================================
 664 /*!
 665  *
 666  */
 667 //=============================================================================
 668 bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 669                                                       Adaptor3d_Curve& theC3d,
 670                                                       double           theLength,
 671                                                       double           theFirstU,
 672                                                       double           theLastU,
 673                                                       list<double> &   theParams,
 674                                                       const bool       theReverse,
 675                                                       bool             theConsiderPropagation)
 676 {
 677   theParams.clear();
 678
 679   double f = theFirstU, l = theLastU;
 680
 681   // Propagation Of Distribution
 682   //
 683   if ( !_mainEdge.IsNull() && _hypType == DISTRIB_PROPAGATION )
 684   {
 685     TopoDS_Edge mainEdge = TopoDS::Edge( _mainEdge ); // should not be a reference!
 686     _gen->Compute( theMesh, mainEdge, SMESH_Gen::SHAPE_ONLY_UPWARD );
 687
 688     SMESHDS_SubMesh* smDS = theMesh.GetMeshDS()->MeshElements( mainEdge );
 689     if ( !smDS )
 690       return error("No mesh on the source edge of Propagation Of Distribution");
 691     if ( smDS->NbNodes() < 1 )
 692       return true; // 1 segment
 693
 694     map< double, const SMDS_MeshNode* > mainEdgeParamsOfNodes;
 695     if ( ! SMESH_Algo::GetSortedNodesOnEdge( theMesh.GetMeshDS(), mainEdge, _quadraticMesh,
 696                                              mainEdgeParamsOfNodes, SMDSAbs_Edge ))
 697       return error("Bad node parameters on the source edge of Propagation Of Distribution");
 698     vector< double > segLen( mainEdgeParamsOfNodes.size() - 1 );
 699     double totalLen = 0;
 700     BRepAdaptor_Curve mainEdgeCurve( mainEdge );
 701     map< double, const SMDS_MeshNode* >::iterator
 702       u_n2 = mainEdgeParamsOfNodes.begin(), u_n1 = u_n2++;
 703     for ( size_t i = 1; i < mainEdgeParamsOfNodes.size(); ++i, ++u_n1, ++u_n2 )
 704     {
 705       segLen[ i-1 ] = GCPnts_AbscissaPoint::Length( mainEdgeCurve,
 706                                                     u_n1->first,
 707                                                     u_n2->first);
 708       totalLen += segLen[ i-1 ];
 709     }
 710     for ( size_t i = 0; i < segLen.size(); ++i )
 711       segLen[ i ] *= theLength / totalLen;
 712
 713     size_t  iSeg = theReverse ? segLen.size()-1 : 0;
 714     size_t  dSeg = theReverse ? -1 : +1;
 715     double param = theFirstU;
 716     size_t nbParams = 0;
 717     for ( int i = 0, nb = segLen.size()-1; i < nb; ++i, iSeg += dSeg )
 718     {
 719       double tol = Min( Precision::Confusion(), 0.01 * segLen[ iSeg ]);
 720       GCPnts_AbscissaPoint Discret( tol, theC3d, segLen[ iSeg ], param );
 721       if ( !Discret.IsDone() ) break;
 722       param = Discret.Parameter();
 723       theParams.push_back( param );
 724       ++nbParams;
 725     }
 726     if ( nbParams != segLen.size()-1 )
 727       return error( SMESH_Comment("Can't divide into ") << segLen.size() << " segments");
 728
 729     compensateError( segLen[ theReverse ? segLen.size()-1 : 0 ],
 730                      segLen[ theReverse ? 0 : segLen.size()-1 ],
 731                      f, l, theLength, theC3d, theParams, true );
 732     return true;
 733   }
 734
 735
 736   switch( _hypType )
 737   {
 738   case LOCAL_LENGTH:
 739   case MAX_LENGTH:
 740   case NB_SEGMENTS:
 741   {
 742     double eltSize = 1;
 743     int nbSegments;
 744     if ( _hypType == MAX_LENGTH )
 745     {
 746       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 747       if (nbseg <= 0)
 748         nbseg = 1; // degenerated edge
 749       eltSize = theLength / nbseg * ( 1. - 1e-9 );
 750       nbSegments = (int) nbseg;
 751     }
 752     else if ( _hypType == LOCAL_LENGTH )
 753     {
 754       // Local Length hypothesis
 755       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 756
 757       // NPAL17873:
 758       bool isFound = false;
 759       if (theConsiderPropagation && !_mainEdge.IsNull()) // propagated from some other edge
 760       {
 761         // Advanced processing to assure equal number of segments in case of Propagation
 762         SMESH_subMesh* sm = theMesh.GetSubMeshContaining(_mainEdge);
 763         if (sm) {
 764           bool computed = sm->IsMeshComputed();
 765           if (!computed) {
 766             if (sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE) {
 767               _gen->Compute( theMesh, _mainEdge, /*anUpward=*/true);
 768               computed = sm->IsMeshComputed();
 769             }
 770           }
 771           if (computed) {
 772             SMESHDS_SubMesh* smds = sm->GetSubMeshDS();
 773             int       nb_segments = smds->NbElements();
 774             if (nbseg - 1 <= nb_segments && nb_segments <= nbseg + 1) {
 775               isFound = true;
 776               nbseg = nb_segments;
 777             }
 778           }
 779         }
 780       }
 781       if (!isFound) // not found by meshed edge in the propagation chain, use precision
 782       {
 783         double aPrecision = _value[ PRECISION_IND ];
 784         double nbseg_prec = ceil((theLength / _value[ BEG_LENGTH_IND ]) - aPrecision);
 785         if (nbseg_prec == (nbseg - 1)) nbseg--;
 786       }
 787
 788       if (nbseg <= 0)
 789         nbseg = 1;                        // degenerated edge
 790       eltSize = theLength / nbseg;
 791       nbSegments = (int) nbseg;
 792     }
 793     else
 794     {
 795       // Number Of Segments hypothesis
 796       nbSegments = _ivalue[ NB_SEGMENTS_IND ];
 797       if ( nbSegments < 1 )  return false;
 798       if ( nbSegments == 1 ) return true;
 799
 800       switch (_ivalue[ DISTR_TYPE_IND ])
 801       {
 802       case StdMeshers_NumberOfSegments::DT_Scale:
 803         {
 804           double scale = _value[ SCALE_FACTOR_IND ];
 805
 806           if (fabs(scale - 1.0) < Precision::Confusion()) {
 807             // special case to avoid division by zero
 808             for (int i = 1; i < nbSegments; i++) {
 809               double param = f + (l - f) * i / nbSegments;
 810               theParams.push_back( param );
 811             }
 812           } else {
 813             // general case of scale distribution
 814             if ( theReverse )
 815               scale = 1.0 / scale;
 816
 817             double  alpha = pow(scale, 1.0 / (nbSegments - 1));
 818             double factor = (l - f) / (1.0 - pow(alpha, nbSegments));
 819
 820             for (int i = 1; i < nbSegments; i++) {
 821               double param = f + factor * (1.0 - pow(alpha, i));
 822               theParams.push_back( param );
 823             }
 824           }
 825           const double lenFactor = theLength/(l-f);
 826           const double minSegLen = Min( theParams.front() - f, l - theParams.back() );
 827           const double       tol = Min( Precision::Confusion(), 0.01 * minSegLen );
 828           list<double>::iterator u = theParams.begin(), uEnd = theParams.end();
 829           for ( ; u != uEnd; ++u )
 830           {
 831             GCPnts_AbscissaPoint Discret( tol, theC3d, ((*u)-f) * lenFactor, f );
 832             if ( Discret.IsDone() )
 833               *u = Discret.Parameter();
 834           }
 835           return true;
 836         }
 837         break;
 838       case StdMeshers_NumberOfSegments::DT_TabFunc:
 839         {
 840           FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
 841           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 842                                     _ivalue[ NB_SEGMENTS_IND ], func,
 843                                     theParams);
 844         }
 845         break;
 846       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 847         {
 848           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
 849           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 850                                     _ivalue[ NB_SEGMENTS_IND ], func,
 851                                     theParams);
 852         }
 853         break;
 854       case StdMeshers_NumberOfSegments::DT_Regular:
 855         eltSize = theLength / nbSegments;
 856         break;
 857       default:
 858         return false;
 859       }
 860     }
 861
 862     double tol = Min( Precision::Confusion(), 0.01 * eltSize );
 863     GCPnts_UniformAbscissa Discret(theC3d, nbSegments + 1, f, l, tol );
 864     if ( !Discret.IsDone() )
 865       return error( "GCPnts_UniformAbscissa failed");
 866     if ( Discret.NbPoints() < nbSegments + 1 )
 867       Discret.Initialize(theC3d, nbSegments + 2, f, l, tol );
 868
 869     int NbPoints = Min( Discret.NbPoints(), nbSegments + 1 );
 870     for ( int i = 2; i < NbPoints; i++ ) // skip 1st and last points
 871     {
 872       double param = Discret.Parameter(i);
 873       theParams.push_back( param );
 874     }
 875     compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams, true ); // for PAL9899
 876     return true;
 877   }
 878
 879
 880   case BEG_END_LENGTH: {
 881
 882     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = theLength
 883
 884     double a1 = _value[ BEG_LENGTH_IND ];
 885     double an = _value[ END_LENGTH_IND ];
 886     double q  = ( theLength - a1 ) / ( theLength - an );
 887     if ( q < theLength/1e6 || 1.01*theLength < a1 + an)
 888       return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
 889                      "for an edge of length "<<theLength);
 890
 891     double      U1 = theReverse ? l : f;
 892     double      Un = theReverse ? f : l;
 893     double   param = U1;
 894     double eltSize = theReverse ? -a1 : a1;
 895     double     tol = Min( Precision::Confusion(), 0.01 * Min( a1, an ));
 896     while ( 1 ) {
 897       // computes a point on a curve <theC3d> at the distance <eltSize>
 898       // from the point of parameter <param>.
 899       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
 900       if ( !Discret.IsDone() ) break;
 901       param = Discret.Parameter();
 902       if ( f < param && param < l )
 903         theParams.push_back( param );
 904       else
 905         break;
 906       eltSize *= q;
 907     }
 908     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 909     if (theReverse) theParams.reverse(); // NPAL18025
 910     return true;
 911   }
 912
 913   case ARITHMETIC_1D:
 914   {
 915     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = theLength
 916
 917     double a1 = _value[ BEG_LENGTH_IND ];
 918     double an = _value[ END_LENGTH_IND ];
 919     if ( 1.01*theLength < a1 + an )
 920       return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
 921                      "for an edge of length "<<theLength);
 922
 923     double q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
 924     int    n = int(fabs(q) > numeric_limits<double>::min() ? ( 1+( an-a1 )/q ) : ( 1+theLength/a1 ));
 925
 926     double      U1 = theReverse ? l : f;
 927     double      Un = theReverse ? f : l;
 928     double   param = U1;
 929     double eltSize = a1;
 930     double     tol = Min( Precision::Confusion(), 0.01 * Min( a1, an ));
 931     if ( theReverse ) {
 932       eltSize = -eltSize;
 933       q = -q;
 934     }
 935     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 936       // computes a point on a curve <theC3d> at the distance <eltSize>
 937       // from the point of parameter <param>.
 938       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
 939       if ( !Discret.IsDone() ) break;
 940       param = Discret.Parameter();
 941       if ( param > f && param < l )
 942         theParams.push_back( param );
 943       else
 944         break;
 945       eltSize += q;
 946     }
 947     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 948     if ( theReverse ) theParams.reverse(); // NPAL18025
 949
 950     return true;
 951   }
 952
 953   case GEOMETRIC_1D:
 954   {
 955     double a1 = _value[ BEG_LENGTH_IND ], an = 0;
 956     double q  = _value[ END_LENGTH_IND ];
 957
 958     double U1 = theReverse ? l : f;
 959     double Un = theReverse ? f : l;
 960     double param = U1;
 961     double eltSize = a1;
 962     if ( theReverse )
 963       eltSize = -eltSize;
 964
 965     int nbParams = 0;
 966     while ( true ) {
 967       // computes a point on a curve <theC3d> at the distance <eltSize>
 968       // from the point of parameter <param>.
 969       double tol = Min( Precision::Confusion(), 0.01 * eltSize );
 970       GCPnts_AbscissaPoint Discret( tol, theC3d, eltSize, param );
 971       if ( !Discret.IsDone() ) break;
 972       param = Discret.Parameter();
 973       if ( f < param && param < l )
 974         theParams.push_back( param );
 975       else
 976         break;
 977       an = eltSize;
 978       eltSize *= q;
 979       ++nbParams;
 980       if ( q < 1. && eltSize < 1e-100 )
 981         return error("Too small common ratio causes too many segments");
 982     }
 983     if ( nbParams > 1 )
 984     {
 985       if ( Abs( param - Un ) < 0.2 * Abs( param - theParams.back() ))
 986       {
 987         compensateError( a1, Abs(eltSize), U1, Un, theLength, theC3d, theParams );
 988       }
 989       else if ( Abs( Un - theParams.back() ) <
 990                 0.2 * Abs( theParams.back() - *(++theParams.rbegin())))
 991       {
 992         theParams.pop_back();
 993         compensateError( a1, Abs(an), U1, Un, theLength, theC3d, theParams );
 994       }
 995     }
 996     if (theReverse) theParams.reverse(); // NPAL18025
 997
 998     return true;
 999   }
1000
1001   case FIXED_POINTS_1D:
1002   {
1003     const std::vector<double>& aPnts = _fpHyp->GetPoints();
1004     std::vector<int>          nbsegs = _fpHyp->GetNbSegments();
1005
1006     // sort normalized params, taking into account theReverse
1007     TColStd_SequenceOfReal Params;
1008     double tol = 1e-7 / theLength; // GCPnts_UniformAbscissa allows u2-u1 > 1e-7
1009     for ( size_t i = 0; i < aPnts.size(); i++ )
1010     {
1011       if( aPnts[i] < tol || aPnts[i] > 1 - tol )
1012         continue;
1013       double u = theReverse ? ( 1 - aPnts[i] ) : aPnts[i];
1014       int    j = 1;
1015       bool IsExist = false;
1016       for ( ; j <= Params.Length(); j++ ) {
1017         if ( Abs( u - Params.Value(j) ) < tol ) {
1018           IsExist = true;
1019           break;
1020         }
1021         if ( u < Params.Value(j) ) break;
1022       }
1023       if ( !IsExist ) Params.InsertBefore( j, u );
1024     }
1025
1026     // transform normalized Params into real ones
1027     std::vector< double > uVec( Params.Length() + 2 );
1028     uVec[ 0 ] = theFirstU;
1029     double abscissa;
1030     for ( int i = 1; i <= Params.Length(); i++ )
1031     {
1032       abscissa = Params( i ) * theLength;
1033       tol      = Min( Precision::Confusion(), 0.01 * abscissa );
1034       GCPnts_AbscissaPoint APnt( tol, theC3d, abscissa, theFirstU );
1035       if ( !APnt.IsDone() )
1036         return error( "GCPnts_AbscissaPoint failed");
1037       uVec[ i ] = APnt.Parameter();
1038     }
1039     uVec.back() = theLastU;
1040
1041     // divide segments
1042     if ( theReverse )
1043     {
1044       if ((int) nbsegs.size() > Params.Length() + 1 )
1045         nbsegs.resize( Params.Length() + 1 );
1046       std::reverse( nbsegs.begin(), nbsegs.end() );
1047     }
1048     if ( nbsegs.empty() )
1049     {
1050       nbsegs.push_back( 1 );
1051     }
1052     Params.InsertBefore( 1, 0.0 );
1053     Params.Append( 1.0 );
1054     double eltSize, segmentSize, par1, par2;
1055     for ( size_t i = 0; i < uVec.size()-1; i++ )
1056     {
1057       par1 = uVec[ i   ];
1058       par2 = uVec[ i+1 ];
1059       int nbseg = ( i < nbsegs.size() ) ? nbsegs[i] : nbsegs[0];
1060       if ( nbseg == 1 )
1061       {
1062         theParams.push_back( par2 );
1063       }
1064       else
1065       {
1066         segmentSize = ( Params( i+2 ) - Params( i+1 )) * theLength;
1067         eltSize     = segmentSize / nbseg;
1068         tol         = Min( Precision::Confusion(), 0.01 * eltSize );
1069         GCPnts_UniformAbscissa Discret( theC3d, eltSize, par1, par2, tol );
1070         if ( !Discret.IsDone() )
1071           return error( "GCPnts_UniformAbscissa failed");
1072         if ( Discret.NbPoints() < nbseg + 1 ) {
1073           eltSize = segmentSize / ( nbseg + 0.5 );
1074           Discret.Initialize( theC3d, eltSize, par1, par2, tol );
1075         }
1076         int NbPoints = Discret.NbPoints();
1077         for ( int i = 2; i <= NbPoints; i++ ) {
1078           double param = Discret.Parameter(i);
1079           theParams.push_back( param );
1080         }
1081       }
1082     }
1083     theParams.pop_back();
1084
1085     return true;
1086   }
1087
1088   case DEFLECTION:
1089   {
1090     GCPnts_UniformDeflection Discret( theC3d, _value[ DEFLECTION_IND ], f, l, true );
1091     if ( !Discret.IsDone() )
1092       return false;
1093
1094     int NbPoints = Discret.NbPoints();
1095     for ( int i = 2; i < NbPoints; i++ )
1096     {
1097       double param = Discret.Parameter(i);
1098       theParams.push_back( param );
1099     }
1100     return true;
1101   }
1102
1103   default:;
1104   }
1105
1106   return false;
1107 }
1108
1109 //=============================================================================
1110 /*!
1111  *
1112  */
1113 //=============================================================================
1114
1115 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & theShape)
1116 {
1117   if ( _hypType == NONE )
1118     return false;
1119
1120   if ( _hypType == ADAPTIVE )
1121   {
1122     _adaptiveHyp->GetAlgo()->InitComputeError();
1123     _adaptiveHyp->GetAlgo()->Compute( theMesh, theShape );
1124     return error( _adaptiveHyp->GetAlgo()->GetComputeError() );
1125   }
1126
1127   SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
1128
1129   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
1130   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
1131   int shapeID = meshDS->ShapeToIndex( E );
1132
1133   double f, l;
1134   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
1135
1136   TopoDS_Vertex VFirst, VLast;
1137   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
1138
1139   ASSERT(!VFirst.IsNull());
1140   ASSERT(!VLast.IsNull());
1141   const SMDS_MeshNode * nFirst = SMESH_Algo::VertexNode( VFirst, meshDS );
1142   const SMDS_MeshNode *  nLast = SMESH_Algo::VertexNode( VLast,  meshDS );
1143   if ( !nFirst || !nLast )
1144     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
1145
1146   // remove elements created by e.g. pattern mapping (PAL21999)
1147   // CLEAN event is incorrectly ptopagated seemingly due to Propagation hyp
1148   // so TEMPORARY solution is to clean the submesh manually
1149   if (SMESHDS_SubMesh * subMeshDS = meshDS->MeshElements(theShape))
1150   {
1151     SMDS_ElemIteratorPtr ite = subMeshDS->GetElements();
1152     while (ite->more())
1153       meshDS->RemoveFreeElement(ite->next(), subMeshDS);
1154     SMDS_NodeIteratorPtr itn = subMeshDS->GetNodes();
1155     while (itn->more()) {
1156       const SMDS_MeshNode * node = itn->next();
1157       if ( node->NbInverseElements() == 0 )
1158         meshDS->RemoveFreeNode(node, subMeshDS);
1159       else
1160         meshDS->RemoveNode(node);
1161     }
1162   }
1163
1164   double length = EdgeLength( E );
1165   if ( !Curve.IsNull() && length > 0 )
1166   {
1167     list< double > params;
1168     bool reversed = false;
1169     if ( theMesh.GetShapeToMesh().ShapeType() >= TopAbs_WIRE && _revEdgesIDs.empty() ) {
1170       // if the shape to mesh is WIRE or EDGE
1171       reversed = ( EE.Orientation() == TopAbs_REVERSED );
1172     }
1173     if ( !_mainEdge.IsNull() ) {
1174       // take into account reversing the edge the hypothesis is propagated from
1175       // (_mainEdge.Orientation() marks mutual orientation of EDGEs in propagation chain)
1176       reversed = ( _mainEdge.Orientation() == TopAbs_REVERSED );
1177       if ( _hypType != DISTRIB_PROPAGATION ) {
1178         int mainID = meshDS->ShapeToIndex(_mainEdge);
1179         if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), mainID) != _revEdgesIDs.end())
1180           reversed = !reversed;
1181       }
1182     }
1183     // take into account this edge reversing
1184     if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), shapeID) != _revEdgesIDs.end())
1185       reversed = !reversed;
1186
1187     BRepAdaptor_Curve C3d( E );
1188     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, reversed, true )) {
1189       return false;
1190     }
1191     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
1192
1193     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
1194     // only internal nodes receive an edge position with param on curve
1195
1196     const SMDS_MeshNode * nPrev = nFirst;
1197     double parPrev = f;
1198     double parLast = l;
1199
1200     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
1201       double param = *itU;
1202       gp_Pnt P = Curve->Value(param);
1203
1204       //Add the Node in the DataStructure
1205       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
1206       meshDS->SetNodeOnEdge(node, shapeID, param);
1207
1208       if(_quadraticMesh) {
1209         // create medium node
1210         double prm = ( parPrev + param )/2;
1211         gp_Pnt  PM = Curve->Value(prm);
1212         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
1213         meshDS->SetNodeOnEdge(NM, shapeID, prm);
1214         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node, NM);
1215         meshDS->SetMeshElementOnShape(edge, shapeID);
1216       }
1217       else {
1218         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node);
1219         meshDS->SetMeshElementOnShape(edge, shapeID);
1220       }
1221
1222       nPrev   = node;
1223       parPrev = param;
1224     }
1225     if(_quadraticMesh) {
1226       double prm = ( parPrev + parLast )/2;
1227       gp_Pnt PM = Curve->Value(prm);
1228       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
1229       meshDS->SetNodeOnEdge(NM, shapeID, prm);
1230       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast, NM);
1231       meshDS->SetMeshElementOnShape(edge, shapeID);
1232     }
1233     else {
1234       SMDS_MeshEdge* edge = meshDS->AddEdge(nPrev, nLast);
1235       meshDS->SetMeshElementOnShape(edge, shapeID);
1236     }
1237   }
1238   else
1239   {
1240     // Edge is a degenerated Edge : We put n = 5 points on the edge.
1241     const int NbPoints = 5;
1242     BRep_Tool::Range( E, f, l ); // PAL15185
1243     double du = (l - f) / (NbPoints - 1);
1244
1245     gp_Pnt P = BRep_Tool::Pnt(VFirst);
1246
1247     const SMDS_MeshNode * nPrev = nFirst;
1248     for (int i = 2; i < NbPoints; i++) {
1249       double param = f + (i - 1) * du;
1250       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
1251       if(_quadraticMesh) {
1252         // create medium node
1253         double prm = param - du/2.;
1254         SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
1255         meshDS->SetNodeOnEdge(NM, shapeID, prm);
1256         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node, NM);
1257         meshDS->SetMeshElementOnShape(edge, shapeID);
1258       }
1259       else {
1260         SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, node);
1261         meshDS->SetMeshElementOnShape(edge, shapeID);
1262       }
1263       meshDS->SetNodeOnEdge(node, shapeID, param);
1264       nPrev = node;
1265     }
1266     if(_quadraticMesh) {
1267       // create medium node
1268       double prm = l - du/2.;
1269       SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
1270       meshDS->SetNodeOnEdge(NM, shapeID, prm);
1271       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast, NM);
1272       meshDS->SetMeshElementOnShape(edge, shapeID);
1273     }
1274     else {
1275       SMDS_MeshEdge * edge = meshDS->AddEdge(nPrev, nLast);
1276       meshDS->SetMeshElementOnShape(edge, shapeID);
1277     }
1278   }
1279   return true;
1280 }
1281
1282
1283 //=============================================================================
1284 /*!
1285  *
1286  */
1287 //=============================================================================
1288
1289 bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh &         theMesh,
1290                                      const TopoDS_Shape & theShape,
1291                                      MapShapeNbElems&     theResMap)
1292 {
1293   if ( _hypType == NONE )
1294     return false;
1295
1296   if ( _hypType == ADAPTIVE )
1297   {
1298     _adaptiveHyp->GetAlgo()->InitComputeError();
1299     _adaptiveHyp->GetAlgo()->Evaluate( theMesh, theShape, theResMap );
1300     return error( _adaptiveHyp->GetAlgo()->GetComputeError() );
1301   }
1302
1303   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
1304   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
1305
1306   double f, l;
1307   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
1308
1309   TopoDS_Vertex VFirst, VLast;
1310   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
1311
1312   ASSERT(!VFirst.IsNull());
1313   ASSERT(!VLast.IsNull());
1314
1315   std::vector<int> aVec(SMDSEntity_Last,0);
1316
1317   double length = EdgeLength( E );
1318   if ( !Curve.IsNull() && length > 0 )
1319   {
1320     list< double > params;
1321     BRepAdaptor_Curve C3d( E );
1322     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, false, true )) {
1323       SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1324       theResMap.insert(std::make_pair(sm,aVec));
1325       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1326       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1327       return false;
1328     }
1329     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
1330
1331     if(_quadraticMesh) {
1332       aVec[SMDSEntity_Node     ] = 2*params.size() + 1;
1333       aVec[SMDSEntity_Quad_Edge] = params.size() + 1;
1334     }
1335     else {
1336       aVec[SMDSEntity_Node] = params.size();
1337       aVec[SMDSEntity_Edge] = params.size() + 1;
1338     }
1339
1340   }
1341   else {
1342     // Edge is a degenerated Edge : We put n = 5 points on the edge.
1343     if ( _quadraticMesh ) {
1344       aVec[SMDSEntity_Node     ] = 11;
1345       aVec[SMDSEntity_Quad_Edge] = 6;
1346     }
1347     else {
1348       aVec[SMDSEntity_Node] = 5;
1349       aVec[SMDSEntity_Edge] = 6;
1350     }
1351   }
1352
1353   SMESH_subMesh * sm = theMesh.GetSubMesh( theShape );
1354   theResMap.insert( std::make_pair( sm, aVec ));
1355
1356   return true;
1357 }
1358
1359
1360 //=============================================================================
1361 /*!
1362  *  See comments in SMESH_Algo.cxx
1363  */
1364 //=============================================================================
1365
1366 const list <const SMESHDS_Hypothesis *> &
1367 StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
1368                                          const TopoDS_Shape & aShape,
1369                                          const bool           ignoreAuxiliary)
1370 {
1371   _usedHypList.clear();
1372   _mainEdge.Nullify();
1373
1374   SMESH_HypoFilter auxiliaryFilter( SMESH_HypoFilter::IsAuxiliary() );
1375   const SMESH_HypoFilter* compatibleFilter = GetCompatibleHypoFilter(/*ignoreAux=*/true );
1376
1377   // get non-auxiliary assigned directly to aShape
1378   int nbHyp = aMesh.GetHypotheses( aShape, *compatibleFilter, _usedHypList, false );
1379
1380   if (nbHyp == 0 && aShape.ShapeType() == TopAbs_EDGE)
1381   {
1382     // Check, if propagated from some other edge
1383     bool isPropagOfDistribution = false;
1384     _mainEdge = StdMeshers_Propagation::GetPropagationSource( aMesh, aShape,
1385                                                               isPropagOfDistribution );
1386     if ( !_mainEdge.IsNull() )
1387     {
1388       if ( isPropagOfDistribution )
1389         _hypType = DISTRIB_PROPAGATION;
1390       // Propagation of 1D hypothesis from <aMainEdge> on this edge;
1391       // get non-auxiliary assigned to _mainEdge
1392       nbHyp = aMesh.GetHypotheses( _mainEdge, *compatibleFilter, _usedHypList, true );
1393     }
1394   }
1395
1396   if (nbHyp == 0) // nothing propagated nor assigned to aShape
1397   {
1398     SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
1399     nbHyp = _usedHypList.size();
1400   }
1401   else
1402   {
1403     // get auxiliary hyps from aShape
1404     aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
1405   }
1406   if ( nbHyp > 1 && ignoreAuxiliary )
1407     _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
1408
1409   return _usedHypList;
1410 }
1411
1412 //================================================================================
1413 /*!
1414  * \brief Pass CancelCompute() to a child algorithm
1415  */
1416 //================================================================================
1417
1418 void StdMeshers_Regular_1D::CancelCompute()
1419 {
1420   SMESH_Algo::CancelCompute();
1421   if ( _hypType == ADAPTIVE )
1422     _adaptiveHyp->GetAlgo()->CancelCompute();
1423 }