src/StdMeshers/StdMeshers_Regular_1D.cxx

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