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         _fpHyp = 0;
  82
  83         _compatibleHypothesis.push_back("LocalLength");
  84         _compatibleHypothesis.push_back("MaxLength");
  85         _compatibleHypothesis.push_back("NumberOfSegments");
  86         _compatibleHypothesis.push_back("StartEndLength");
  87         _compatibleHypothesis.push_back("Deflection1D");
  88         _compatibleHypothesis.push_back("Arithmetic1D");
  89         _compatibleHypothesis.push_back("FixedPoints1D");
  90         _compatibleHypothesis.push_back("AutomaticLength");
  91
  92         _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
  93         _compatibleHypothesis.push_back("Propagation"); // auxiliary !!!
  94 }
  95
  96 //=============================================================================
  97 /*!
  98  *
  99  */
 100 //=============================================================================
 101
 102 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
 103 {
 104 }
 105
 106 //=============================================================================
 107 /*!
 108  *
 109  */
 110 //=============================================================================
 111
 112 bool StdMeshers_Regular_1D::CheckHypothesis
 113                          (SMESH_Mesh&                          aMesh,
 114                           const TopoDS_Shape&                  aShape,
 115                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 116 {
 117   _hypType = NONE;
 118   _quadraticMesh = false;
 119
 120   const bool ignoreAuxiliaryHyps = false;
 121   const list <const SMESHDS_Hypothesis * > & hyps =
 122     GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
 123
 124   // find non-auxiliary hypothesis
 125   const SMESHDS_Hypothesis *theHyp = 0;
 126   list <const SMESHDS_Hypothesis * >::const_iterator h = hyps.begin();
 127   for ( ; h != hyps.end(); ++h ) {
 128     if ( static_cast<const SMESH_Hypothesis*>(*h)->IsAuxiliary() ) {
 129       if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 )
 130         _quadraticMesh = true;
 131     }
 132     else {
 133       if ( !theHyp )
 134         theHyp = *h; // use only the first non-auxiliary hypothesis
 135     }
 136   }
 137
 138   if ( !theHyp )
 139   {
 140     aStatus = SMESH_Hypothesis::HYP_MISSING;
 141     return false;  // can't work without a hypothesis
 142   }
 143
 144   string hypName = theHyp->GetName();
 145
 146   if (hypName == "LocalLength")
 147   {
 148     const StdMeshers_LocalLength * hyp =
 149       dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
 150     ASSERT(hyp);
 151     _value[ BEG_LENGTH_IND ] = hyp->GetLength();
 152     _value[ PRECISION_IND ] = hyp->GetPrecision();
 153     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 154     _hypType = LOCAL_LENGTH;
 155     aStatus = SMESH_Hypothesis::HYP_OK;
 156   }
 157
 158   else if (hypName == "MaxLength")
 159   {
 160     const StdMeshers_MaxLength * hyp =
 161       dynamic_cast <const StdMeshers_MaxLength * >(theHyp);
 162     ASSERT(hyp);
 163     _value[ BEG_LENGTH_IND ] = hyp->GetLength();
 164     if ( hyp->GetUsePreestimatedLength() ) {
 165       if ( int nbSeg = aMesh.GetGen()->GetBoundaryBoxSegmentation() )
 166         _value[ BEG_LENGTH_IND ] = aMesh.GetShapeDiagonalSize() / nbSeg;
 167     }
 168     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 169     _hypType = MAX_LENGTH;
 170     aStatus = SMESH_Hypothesis::HYP_OK;
 171   }
 172
 173   else if (hypName == "NumberOfSegments")
 174   {
 175     const StdMeshers_NumberOfSegments * hyp =
 176       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
 177     ASSERT(hyp);
 178     _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
 179     ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
 180     _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
 181     switch (_ivalue[ DISTR_TYPE_IND ])
 182     {
 183     case StdMeshers_NumberOfSegments::DT_Scale:
 184       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
 185       _revEdgesIDs = hyp->GetReversedEdges();
 186       break;
 187     case StdMeshers_NumberOfSegments::DT_TabFunc:
 188       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
 189       _revEdgesIDs = hyp->GetReversedEdges();
 190       break;
 191     case StdMeshers_NumberOfSegments::DT_ExprFunc:
 192       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
 193       _revEdgesIDs = hyp->GetReversedEdges();
 194       break;
 195     case StdMeshers_NumberOfSegments::DT_Regular:
 196       break;
 197     default:
 198       ASSERT(0);
 199       break;
 200     }
 201     if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
 202         _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
 203         _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
 204     _hypType = NB_SEGMENTS;
 205     aStatus = SMESH_Hypothesis::HYP_OK;
 206   }
 207
 208   else if (hypName == "Arithmetic1D")
 209   {
 210     const StdMeshers_Arithmetic1D * hyp =
 211       dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
 212     ASSERT(hyp);
 213     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 214     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 215     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 216     _hypType = ARITHMETIC_1D;
 217
 218     _revEdgesIDs = hyp->GetReversedEdges();
 219
 220     aStatus = SMESH_Hypothesis::HYP_OK;
 221   }
 222
 223   else if (hypName == "FixedPoints1D") {
 224     _fpHyp = dynamic_cast <const StdMeshers_FixedPoints1D*>(theHyp);
 225     ASSERT(_fpHyp);
 226     _hypType = FIXED_POINTS_1D;
 227
 228     _revEdgesIDs = _fpHyp->GetReversedEdges();
 229
 230     aStatus = SMESH_Hypothesis::HYP_OK;
 231   }
 232
 233   else if (hypName == "StartEndLength")
 234   {
 235     const StdMeshers_StartEndLength * hyp =
 236       dynamic_cast <const StdMeshers_StartEndLength * >(theHyp);
 237     ASSERT(hyp);
 238     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 239     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 240     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 241     _hypType = BEG_END_LENGTH;
 242
 243     _revEdgesIDs = hyp->GetReversedEdges();
 244
 245     aStatus = SMESH_Hypothesis::HYP_OK;
 246   }
 247
 248   else if (hypName == "Deflection1D")
 249   {
 250     const StdMeshers_Deflection1D * hyp =
 251       dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
 252     ASSERT(hyp);
 253     _value[ DEFLECTION_IND ] = hyp->GetDeflection();
 254     ASSERT( _value[ DEFLECTION_IND ] > 0 );
 255     _hypType = DEFLECTION;
 256     aStatus = SMESH_Hypothesis::HYP_OK;
 257   }
 258
 259   else if (hypName == "AutomaticLength")
 260   {
 261     StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
 262       (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
 263     ASSERT(hyp);
 264     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 265 //     _value[ BEG_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 266 //     _value[ END_LENGTH_IND ] = Precision::Confusion(); // ?? or set to zero?
 267     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 268     _hypType = MAX_LENGTH;
 269     aStatus = SMESH_Hypothesis::HYP_OK;
 270   }
 271   else
 272     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 273
 274   return ( _hypType != NONE );
 275 }
 276
 277 static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
 278                                double length, bool theReverse,
 279                                int nbSeg, Function& func,
 280                                list<double>& theParams)
 281 {
 282   // never do this way
 283   //OSD::SetSignal( true );
 284
 285   if (nbSeg <= 0)
 286     return false;
 287
 288   MESSAGE( "computeParamByFunc" );
 289
 290   int nbPnt = 1 + nbSeg;
 291   vector<double> x(nbPnt, 0.);
 292
 293   if (!buildDistribution(func, 0.0, 1.0, nbSeg, x, 1E-4))
 294      return false;
 295
 296   MESSAGE( "Points:\n" );
 297   char buf[1024];
 298   for ( int i=0; i<=nbSeg; i++ )
 299   {
 300     sprintf(  buf, "%f\n", float(x[i] ) );
 301     MESSAGE( buf );
 302   }
 303
 304
 305
 306   // apply parameters in range [0,1] to the space of the curve
 307   double prevU = first;
 308   double sign = 1.;
 309   if (theReverse)
 310   {
 311     prevU = last;
 312     sign = -1.;
 313   }
 314   for( int i = 1; i < nbSeg; i++ )
 315   {
 316     double curvLength = length * (x[i] - x[i-1]) * sign;
 317     GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
 318     if ( !Discret.IsDone() )
 319       return false;
 320     double U = Discret.Parameter();
 321     if ( U > first && U < last )
 322       theParams.push_back( U );
 323     else
 324       return false;
 325     prevU = U;
 326   }
 327   return true;
 328 }
 329
 330
 331 //================================================================================
 332 /*!
 333  * \brief adjust internal node parameters so that the last segment length == an
 334   * \param a1 - the first segment length
 335   * \param an - the last segment length
 336   * \param U1 - the first edge parameter
 337   * \param Un - the last edge parameter
 338   * \param length - the edge length
 339   * \param C3d - the edge curve
 340   * \param theParams - internal node parameters to adjust
 341   * \param adjustNeighbors2an - to adjust length of segments next to the last one
 342   *  and not to remove parameters
 343  */
 344 //================================================================================
 345
 346 static void compensateError(double a1, double an,
 347                             double U1, double Un,
 348                             double            length,
 349                             Adaptor3d_Curve&  C3d,
 350                             list<double> &    theParams,
 351                             bool              adjustNeighbors2an = false)
 352 {
 353   int i, nPar = theParams.size();
 354   if ( a1 + an < length && nPar > 1 )
 355   {
 356     bool reverse = ( U1 > Un );
 357     GCPnts_AbscissaPoint Discret(C3d, reverse ? an : -an, Un);
 358     if ( !Discret.IsDone() )
 359       return;
 360     double Utgt = Discret.Parameter(); // target value of the last parameter
 361     list<double>::reverse_iterator itU = theParams.rbegin();
 362     double Ul = *itU++; // real value of the last parameter
 363     double dUn = Utgt - Ul; // parametric error of <an>
 364     if ( Abs(dUn) <= Precision::Confusion() )
 365       return;
 366     double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
 367     if ( adjustNeighbors2an || Abs(dUn) < 0.5 * dU ) { // last segment is a bit shorter than it should
 368       // move the last parameter to the edge beginning
 369     }
 370     else {  // last segment is much shorter than it should -> remove the last param and
 371       theParams.pop_back(); nPar--; // move the rest points toward the edge end
 372       dUn = Utgt - theParams.back();
 373     }
 374
 375     double q  = dUn / ( nPar - 1 );
 376     if ( !adjustNeighbors2an ) {
 377       for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
 378         (*itU) += dUn;
 379         dUn -= q;
 380       }
 381     }
 382     else {
 383       theParams.back() += dUn;
 384       double sign = reverse ? -1 : 1;
 385       double prevU = theParams.back();
 386       itU = theParams.rbegin();
 387       for ( ++itU, i = 2; i < nPar; ++itU, i++ ) {
 388         double newU = *itU + dUn;
 389         if ( newU*sign < prevU*sign ) {
 390           prevU = *itU = newU;
 391           dUn -= q;
 392         }
 393         else { // set U between prevU and next valid param
 394           list<double>::reverse_iterator itU2 = itU;
 395           ++itU2;
 396           int nb = 2;
 397           while ( (*itU2)*sign > prevU*sign ) {
 398             ++itU2; ++nb;
 399           }
 400           dU = ( *itU2 - prevU ) / nb;
 401           while ( itU != itU2 ) {
 402             *itU += dU; ++itU;
 403           }
 404           break;
 405         }
 406       }
 407     }
 408   }
 409 }
 410
 411 //================================================================================
 412 /*!
 413  * \brief Class used to clean mesh on edges when 0D hyp modified.
 414  * Common approach doesn't work when 0D algo is missing because the 0D hyp is
 415  * considered as not participating in computation whereas it is used by 1D algo.
 416  */
 417 //================================================================================
 418
 419 // struct VertexEventListener : public SMESH_subMeshEventListener
 420 // {
 421 //   VertexEventListener():SMESH_subMeshEventListener(0) // won't be deleted by submesh
 422 //   {}
 423 //   /*!
 424 //    * \brief Clean mesh on edges
 425 //    * \param event - algo_event or compute_event itself (of SMESH_subMesh)
 426 //    * \param eventType - ALGO_EVENT or COMPUTE_EVENT (of SMESH_subMesh)
 427 //    * \param subMesh - the submesh where the event occures
 428 //    */
 429 //   void ProcessEvent(const int event, const int eventType, SMESH_subMesh* subMesh,
 430 //                     EventListenerData*, const SMESH_Hypothesis*)
 431 //   {
 432 //     if ( eventType == SMESH_subMesh::ALGO_EVENT) // all algo events
 433 //     {
 434 //       subMesh->ComputeStateEngine( SMESH_subMesh::MODIF_ALGO_STATE );
 435 //     }
 436 //   }
 437 // }; // struct VertexEventListener
 438
 439 //=============================================================================
 440 /*!
 441  * \brief Sets event listener to vertex submeshes
 442  * \param subMesh - submesh where algo is set
 443  *
 444  * This method is called when a submesh gets HYP_OK algo_state.
 445  * After being set, event listener is notified on each event of a submesh.
 446  */
 447 //=============================================================================
 448
 449 void StdMeshers_Regular_1D::SetEventListener(SMESH_subMesh* subMesh)
 450 {
 451   StdMeshers_Propagation::SetPropagationMgr( subMesh );
 452 }
 453
 454 //=============================================================================
 455 /*!
 456  * \brief Do nothing
 457  * \param subMesh - restored submesh
 458  *
 459  * This method is called only if a submesh has HYP_OK algo_state.
 460  */
 461 //=============================================================================
 462
 463 void StdMeshers_Regular_1D::SubmeshRestored(SMESH_subMesh* subMesh)
 464 {
 465 }
 466
 467 //=============================================================================
 468 /*!
 469  * \brief Return StdMeshers_SegmentLengthAroundVertex assigned to vertex
 470  */
 471 //=============================================================================
 472
 473 const StdMeshers_SegmentLengthAroundVertex*
 474 StdMeshers_Regular_1D::getVertexHyp(SMESH_Mesh &          theMesh,
 475                                     const TopoDS_Vertex & theV)
 476 {
 477   static SMESH_HypoFilter filter( SMESH_HypoFilter::HasName("SegmentAroundVertex_0D"));
 478   if ( const SMESH_Hypothesis * h = theMesh.GetHypothesis( theV, filter, true ))
 479   {
 480     SMESH_Algo* algo = const_cast< SMESH_Algo* >( static_cast< const SMESH_Algo* > ( h ));
 481     const list <const SMESHDS_Hypothesis *> & hypList = algo->GetUsedHypothesis( theMesh, theV, 0 );
 482     if ( !hypList.empty() && string("SegmentLengthAroundVertex") == hypList.front()->GetName() )
 483       return static_cast<const StdMeshers_SegmentLengthAroundVertex*>( hypList.front() );
 484   }
 485   return 0;
 486 }
 487
 488 //================================================================================
 489 /*!
 490  * \brief Tune parameters to fit "SegmentLengthAroundVertex" hypothesis
 491   * \param theC3d - wire curve
 492   * \param theLength - curve length
 493   * \param theParameters - internal nodes parameters to modify
 494   * \param theVf - 1st vertex
 495   * \param theVl - 2nd vertex
 496  */
 497 //================================================================================
 498
 499 void StdMeshers_Regular_1D::redistributeNearVertices (SMESH_Mesh &          theMesh,
 500                                                       Adaptor3d_Curve &     theC3d,
 501                                                       double                theLength,
 502                                                       std::list< double > & theParameters,
 503                                                       const TopoDS_Vertex & theVf,
 504                                                       const TopoDS_Vertex & theVl)
 505 {
 506   double f = theC3d.FirstParameter(), l = theC3d.LastParameter();
 507   int nPar = theParameters.size();
 508   for ( int isEnd1 = 0; isEnd1 < 2; ++isEnd1 )
 509   {
 510     const TopoDS_Vertex & V = isEnd1 ? theVf : theVl;
 511     const StdMeshers_SegmentLengthAroundVertex* hyp = getVertexHyp (theMesh, V );
 512     if ( hyp ) {
 513       double vertexLength = hyp->GetLength();
 514       if ( vertexLength > theLength / 2.0 )
 515         continue;
 516       if ( isEnd1 ) { // to have a segment of interest at end of theParameters
 517         theParameters.reverse();
 518         std::swap( f, l );
 519       }
 520       if ( _hypType == NB_SEGMENTS )
 521       {
 522         compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 523       }
 524       else if ( nPar <= 3 )
 525       {
 526         if ( !isEnd1 )
 527           vertexLength = -vertexLength;
 528         GCPnts_AbscissaPoint Discret(theC3d, vertexLength, l);
 529         if ( Discret.IsDone() ) {
 530           if ( nPar == 0 )
 531             theParameters.push_back( Discret.Parameter());
 532           else {
 533             double L = GCPnts_AbscissaPoint::Length( theC3d, theParameters.back(), l);
 534             if ( vertexLength < L / 2.0 )
 535               theParameters.push_back( Discret.Parameter());
 536             else
 537               compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 538           }
 539         }
 540       }
 541       else
 542       {
 543         // recompute params between the last segment and a middle one.
 544         // find size of a middle segment
 545         int nHalf = ( nPar-1 ) / 2;
 546         list< double >::reverse_iterator itU = theParameters.rbegin();
 547         std::advance( itU, nHalf );
 548         double Um = *itU++;
 549         double Lm = GCPnts_AbscissaPoint::Length( theC3d, Um, *itU);
 550         double L = GCPnts_AbscissaPoint::Length( theC3d, *itU, l);
 551         StdMeshers_Regular_1D algo( *this );
 552         algo._hypType = BEG_END_LENGTH;
 553         algo._value[ BEG_LENGTH_IND ] = Lm;
 554         algo._value[ END_LENGTH_IND ] = vertexLength;
 555         double from = *itU, to = l;
 556         if ( isEnd1 ) {
 557           std::swap( from, to );
 558           std::swap( algo._value[ BEG_LENGTH_IND ], algo._value[ END_LENGTH_IND ]);
 559         }
 560         list<double> params;
 561         if ( algo.computeInternalParameters( theMesh, theC3d, L, from, to, params, false ))
 562         {
 563           if ( isEnd1 ) params.reverse();
 564           while ( 1 + nHalf-- )
 565             theParameters.pop_back();
 566           theParameters.splice( theParameters.end(), params );
 567         }
 568         else
 569         {
 570           compensateError(0, vertexLength, f, l, theLength, theC3d, theParameters, true );
 571         }
 572       }
 573       if ( isEnd1 )
 574         theParameters.reverse();
 575     }
 576   }
 577 }
 578
 579 //=============================================================================
 580 /*!
 581  *
 582  */
 583 //=============================================================================
 584 bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 585                                                       Adaptor3d_Curve& theC3d,
 586                                                       double           theLength,
 587                                                       double           theFirstU,
 588                                                       double           theLastU,
 589                                                       list<double> &   theParams,
 590                                                       const bool       theReverse,
 591                                                       bool             theConsiderPropagation)
 592 {
 593   theParams.clear();
 594
 595   double f = theFirstU, l = theLastU;
 596
 597   switch( _hypType )
 598   {
 599   case LOCAL_LENGTH:
 600   case MAX_LENGTH:
 601   case NB_SEGMENTS: {
 602
 603     double eltSize = 1;
 604     if ( _hypType == MAX_LENGTH )
 605     {
 606       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 607       if (nbseg <= 0)
 608         nbseg = 1;                        // degenerated edge
 609       eltSize = theLength / nbseg;
 610     }
 611     else if ( _hypType == LOCAL_LENGTH )
 612     {
 613       // Local Length hypothesis
 614       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
 615
 616       // NPAL17873:
 617       bool isFound = false;
 618       if (theConsiderPropagation && !_mainEdge.IsNull()) // propagated from some other edge
 619       {
 620         // Advanced processing to assure equal number of segments in case of Propagation
 621         SMESH_subMesh* sm = theMesh.GetSubMeshContaining(_mainEdge);
 622         if (sm) {
 623           bool computed = sm->IsMeshComputed();
 624           if (!computed) {
 625             if (sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE) {
 626               sm->ComputeStateEngine(SMESH_subMesh::COMPUTE);
 627               computed = sm->IsMeshComputed();
 628             }
 629           }
 630           if (computed) {
 631             SMESHDS_SubMesh* smds = sm->GetSubMeshDS();
 632             int nb_segments = smds->NbElements();
 633             if (nbseg - 1 <= nb_segments && nb_segments <= nbseg + 1) {
 634               isFound = true;
 635               nbseg = nb_segments;
 636             }
 637           }
 638         }
 639       }
 640       if (!isFound) // not found by meshed edge in the propagation chain, use precision
 641       {
 642         double aPrecision = _value[ PRECISION_IND ];
 643         double nbseg_prec = ceil((theLength / _value[ BEG_LENGTH_IND ]) - aPrecision);
 644         if (nbseg_prec == (nbseg - 1)) nbseg--;
 645       }
 646
 647       if (nbseg <= 0)
 648         nbseg = 1;                        // degenerated edge
 649       eltSize = theLength / nbseg;
 650     }
 651     else
 652     {
 653       // Number Of Segments hypothesis
 654       int NbSegm = _ivalue[ NB_SEGMENTS_IND ];
 655       if ( NbSegm < 1 )  return false;
 656       if ( NbSegm == 1 ) return true;
 657
 658       switch (_ivalue[ DISTR_TYPE_IND ])
 659       {
 660       case StdMeshers_NumberOfSegments::DT_Scale:
 661         {
 662           double scale = _value[ SCALE_FACTOR_IND ];
 663
 664           if (fabs(scale - 1.0) < Precision::Confusion()) {
 665             // special case to avoid division by zero
 666             for (int i = 1; i < NbSegm; i++) {
 667               double param = f + (l - f) * i / NbSegm;
 668               theParams.push_back( param );
 669             }
 670           } else {
 671             // general case of scale distribution
 672             if ( theReverse )
 673               scale = 1.0 / scale;
 674
 675             double alpha = pow(scale, 1.0 / (NbSegm - 1));
 676             double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
 677
 678             for (int i = 1; i < NbSegm; i++) {
 679               double param = f + factor * (1.0 - pow(alpha, i));
 680               theParams.push_back( param );
 681             }
 682           }
 683           return true;
 684         }
 685         break;
 686       case StdMeshers_NumberOfSegments::DT_TabFunc:
 687         {
 688           FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
 689           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 690                                     _ivalue[ NB_SEGMENTS_IND ], func,
 691                                     theParams);
 692         }
 693         break;
 694       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 695         {
 696           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
 697           return computeParamByFunc(theC3d, f, l, theLength, theReverse,
 698                                     _ivalue[ NB_SEGMENTS_IND ], func,
 699                                     theParams);
 700         }
 701         break;
 702       case StdMeshers_NumberOfSegments::DT_Regular:
 703         eltSize = theLength / _ivalue[ NB_SEGMENTS_IND ];
 704         break;
 705       default:
 706         return false;
 707       }
 708     }
 709     GCPnts_UniformAbscissa Discret(theC3d, eltSize, f, l);
 710     if ( !Discret.IsDone() )
 711       return error( "GCPnts_UniformAbscissa failed");
 712
 713     int NbPoints = Discret.NbPoints();
 714     for ( int i = 2; i < NbPoints; i++ )
 715     {
 716       double param = Discret.Parameter(i);
 717       theParams.push_back( param );
 718     }
 719     compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams ); // for PAL9899
 720     return true;
 721   }
 722
 723   case BEG_END_LENGTH: {
 724
 725     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = theLength
 726
 727     double a1 = _value[ BEG_LENGTH_IND ];
 728     double an = _value[ END_LENGTH_IND ];
 729     double q  = ( theLength - a1 ) / ( theLength - an );
 730
 731     double U1 = theReverse ? l : f;
 732     double Un = theReverse ? f : l;
 733     double param = U1;
 734     double eltSize = theReverse ? -a1 : a1;
 735     while ( 1 ) {
 736       // computes a point on a curve <theC3d> at the distance <eltSize>
 737       // from the point of parameter <param>.
 738       GCPnts_AbscissaPoint Discret( theC3d, eltSize, param );
 739       if ( !Discret.IsDone() ) break;
 740       param = Discret.Parameter();
 741       if ( f < param && param < l )
 742         theParams.push_back( param );
 743       else
 744         break;
 745       eltSize *= q;
 746     }
 747     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 748     if (theReverse) theParams.reverse(); // NPAL18025
 749     return true;
 750   }
 751
 752   case ARITHMETIC_1D: {
 753
 754     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = theLength
 755
 756     double a1 = _value[ BEG_LENGTH_IND ];
 757     double an = _value[ END_LENGTH_IND ];
 758
 759     double  q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
 760     int     n = int( 1 + ( an - a1 ) / q );
 761
 762     double U1 = theReverse ? l : f;
 763     double Un = theReverse ? f : l;
 764     double param = U1;
 765     double eltSize = a1;
 766     if ( theReverse ) {
 767       eltSize = -eltSize;
 768       q = -q;
 769     }
 770     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 771       // computes a point on a curve <theC3d> at the distance <eltSize>
 772       // from the point of parameter <param>.
 773       GCPnts_AbscissaPoint Discret( theC3d, eltSize, param );
 774       if ( !Discret.IsDone() ) break;
 775       param = Discret.Parameter();
 776       if ( param > f && param < l )
 777         theParams.push_back( param );
 778       else
 779         break;
 780       eltSize += q;
 781     }
 782     compensateError( a1, an, U1, Un, theLength, theC3d, theParams );
 783     if (theReverse) theParams.reverse(); // NPAL18025
 784
 785     return true;
 786   }
 787
 788   case FIXED_POINTS_1D: {
 789     const std::vector<double>& aPnts = _fpHyp->GetPoints();
 790     const std::vector<int>& nbsegs = _fpHyp->GetNbSegments();
 791     int i = 0;
 792     TColStd_SequenceOfReal Params;
 793     for(; i<aPnts.size(); i++) {
 794       if( aPnts[i]<0.0001 || aPnts[i]>0.9999 ) continue;
 795       int j=1;
 796       bool IsExist = false;
 797       for(; j<=Params.Length(); j++) {
 798         if( fabs(aPnts[i]-Params.Value(j)) < 1e-4 ) {
 799           IsExist = true;
 800           break;
 801         }
 802         if( aPnts[i]<Params.Value(j) ) break;
 803       }
 804       if(!IsExist) Params.InsertBefore(j,aPnts[i]);
 805     }
 806     double pf, pl, par2, par1, psize;
 807     if (theReverse) {
 808       pf = l;
 809       pl = f;
 810     }
 811     else {
 812       pf = f;
 813       pl = l;
 814     }
 815     psize = pl - pf;
 816     par1 = pf;
 817     //cout<<"aPnts.size() = "<<aPnts.size()<<"  Params.Length() = "
 818     //    <<Params.Length()<<"   nbsegs.size() = "<<nbsegs.size()<<endl;
 819     for(i=0; i<Params.Length(); i++) {
 820       par2 = pf + Params.Value(i+1)*psize;
 821       int nbseg = ( i > nbsegs.size()-1 ) ? nbsegs[0] : nbsegs[i];
 822       double dp = (par2-par1)/nbseg;
 823       int j = 1;
 824       for(; j<=nbseg; j++) {
 825         double param = par1 + dp*j;
 826         theParams.push_back( param );
 827       }
 828       par1 = par2;
 829     }
 830     // add for last
 831     int nbseg = ( nbsegs.size() > Params.Length() ) ? nbsegs[Params.Length()] : nbsegs[0];
 832     double dp = (pl-par1)/nbseg;
 833     int j = 1;
 834     for(; j<nbseg; j++) {
 835       double param = par1 + dp*j;
 836       theParams.push_back( param );
 837     }
 838     if (theReverse) {
 839       theParams.reverse(); // NPAL18025
 840     }
 841     return true;
 842   }
 843
 844   case DEFLECTION: {
 845
 846     GCPnts_UniformDeflection Discret(theC3d, _value[ DEFLECTION_IND ], f, l, true);
 847     if ( !Discret.IsDone() )
 848       return false;
 849
 850     int NbPoints = Discret.NbPoints();
 851     for ( int i = 2; i < NbPoints; i++ )
 852     {
 853       double param = Discret.Parameter(i);
 854       theParams.push_back( param );
 855     }
 856     return true;
 857   }
 858
 859   default:;
 860   }
 861
 862   return false;
 863 }
 864
 865 //=============================================================================
 866 /*!
 867  *
 868  */
 869 //=============================================================================
 870
 871 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & theShape)
 872 {
 873   if ( _hypType == NONE )
 874     return false;
 875
 876   SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
 877
 878   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
 879   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
 880   int shapeID = meshDS->ShapeToIndex( E );
 881
 882   double f, l;
 883   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
 884
 885   TopoDS_Vertex VFirst, VLast;
 886   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
 887
 888   ASSERT(!VFirst.IsNull());
 889   ASSERT(!VLast.IsNull());
 890   const SMDS_MeshNode * idFirst = SMESH_Algo::VertexNode( VFirst, meshDS );
 891   const SMDS_MeshNode * idLast = SMESH_Algo::VertexNode( VLast, meshDS );
 892   if (!idFirst || !idLast)
 893     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
 894
 895   if (!Curve.IsNull())
 896   {
 897     list< double > params;
 898     bool reversed = false;
 899     if ( theMesh.GetShapeToMesh().ShapeType() >= TopAbs_WIRE ) {
 900       reversed = ( EE.Orientation() == TopAbs_REVERSED );
 901     }
 902     if ( !_mainEdge.IsNull() ) {
 903       reversed = ( _mainEdge.Orientation() == TopAbs_REVERSED );
 904     }
 905     else if ( _revEdgesIDs.size() > 0 ) {
 906       for ( int i = 0; i < _revEdgesIDs.size(); i++) {
 907         if ( _revEdgesIDs[i] == shapeID ) {
 908           reversed = !reversed;
 909         }
 910       }
 911     }
 912
 913     BRepAdaptor_Curve C3d( E );
 914     double length = EdgeLength( E );
 915     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, reversed, true )) {
 916       return false;
 917     }
 918     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
 919
 920     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
 921     // only internal nodes receive an edge position with param on curve
 922
 923     const SMDS_MeshNode * idPrev = idFirst;
 924     double parPrev = f;
 925     double parLast = l;
 926
 927     /* NPAL18025
 928     if (reversed) {
 929       idPrev = idLast;
 930       idLast = idFirst;
 931       idFirst = idPrev;
 932       parPrev = l;
 933       parLast = f;
 934     }
 935     */
 936
 937     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
 938       double param = *itU;
 939       gp_Pnt P = Curve->Value(param);
 940
 941       //Add the Node in the DataStructure
 942       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 943       meshDS->SetNodeOnEdge(node, shapeID, param);
 944
 945       if(_quadraticMesh) {
 946         // create medium node
 947         double prm = ( parPrev + param )/2;
 948         gp_Pnt PM = Curve->Value(prm);
 949         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 950         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 951         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 952         meshDS->SetMeshElementOnShape(edge, shapeID);
 953       }
 954       else {
 955         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 956         meshDS->SetMeshElementOnShape(edge, shapeID);
 957       }
 958
 959       idPrev = node;
 960       parPrev = param;
 961     }
 962     if(_quadraticMesh) {
 963       double prm = ( parPrev + parLast )/2;
 964       gp_Pnt PM = Curve->Value(prm);
 965       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 966       meshDS->SetNodeOnEdge(NM, shapeID, prm);
 967       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
 968       meshDS->SetMeshElementOnShape(edge, shapeID);
 969     }
 970     else {
 971       SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
 972       meshDS->SetMeshElementOnShape(edge, shapeID);
 973     }
 974   }
 975   else
 976   {
 977     //MESSAGE("************* Degenerated edge! *****************");
 978
 979     // Edge is a degenerated Edge : We put n = 5 points on the edge.
 980     const int NbPoints = 5;
 981     BRep_Tool::Range( E, f, l ); // PAL15185
 982     double du = (l - f) / (NbPoints - 1);
 983
 984     gp_Pnt P = BRep_Tool::Pnt(VFirst);
 985
 986     const SMDS_MeshNode * idPrev = idFirst;
 987     for (int i = 2; i < NbPoints; i++) {
 988       double param = f + (i - 1) * du;
 989       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 990       if(_quadraticMesh) {
 991         // create medium node
 992         double prm = param - du/2.;
 993         SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
 994         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 995         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 996         meshDS->SetMeshElementOnShape(edge, shapeID);
 997       }
 998       else {
 999         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
1000         meshDS->SetMeshElementOnShape(edge, shapeID);
1001       }
1002       meshDS->SetNodeOnEdge(node, shapeID, param);
1003       idPrev = node;
1004     }
1005     if(_quadraticMesh) {
1006       // create medium node
1007       double prm = l - du/2.;
1008       SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z());
1009       meshDS->SetNodeOnEdge(NM, shapeID, prm);
1010       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
1011       meshDS->SetMeshElementOnShape(edge, shapeID);
1012     }
1013     else {
1014       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
1015       meshDS->SetMeshElementOnShape(edge, shapeID);
1016     }
1017   }
1018   return true;
1019 }
1020
1021
1022 //=============================================================================
1023 /*!
1024  *
1025  */
1026 //=============================================================================
1027
1028 bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh & theMesh,
1029                                      const TopoDS_Shape & theShape,
1030                                      MapShapeNbElems& aResMap)
1031 {
1032   if ( _hypType == NONE )
1033     return false;
1034
1035   SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
1036
1037   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
1038   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
1039   //  int shapeID = meshDS->ShapeToIndex( E );
1040
1041   double f, l;
1042   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
1043
1044   TopoDS_Vertex VFirst, VLast;
1045   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
1046
1047   ASSERT(!VFirst.IsNull());
1048   ASSERT(!VLast.IsNull());
1049
1050   std::vector<int> aVec(SMDSEntity_Last);
1051   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
1052
1053   if (!Curve.IsNull()) {
1054     list< double > params;
1055
1056     BRepAdaptor_Curve C3d( E );
1057     double length = EdgeLength( E );
1058     if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, false, true )) {
1059       SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1060       aResMap.insert(std::make_pair(sm,aVec));
1061       SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1062       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
1063       return false;
1064     }
1065     redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
1066
1067     if(_quadraticMesh) {
1068       aVec[SMDSEntity_Node] = 2*params.size() + 1;
1069       aVec[SMDSEntity_Quad_Edge] = params.size() + 1;
1070     }
1071     else {
1072       aVec[SMDSEntity_Node] = params.size();
1073       aVec[SMDSEntity_Edge] = params.size() + 1;
1074     }
1075
1076   }
1077   else {
1078     //MESSAGE("************* Degenerated edge! *****************");
1079     // Edge is a degenerated Edge : We put n = 5 points on the edge.
1080     if(_quadraticMesh) {
1081       aVec[SMDSEntity_Node] = 11;
1082       aVec[SMDSEntity_Quad_Edge] = 6;
1083     }
1084     else {
1085       aVec[SMDSEntity_Node] = 5;
1086       aVec[SMDSEntity_Edge] = 6;
1087     }
1088   }
1089
1090   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1091   aResMap.insert(std::make_pair(sm,aVec));
1092
1093   return true;
1094 }
1095
1096
1097 //=============================================================================
1098 /*!
1099  *  See comments in SMESH_Algo.cxx
1100  */
1101 //=============================================================================
1102
1103 const list <const SMESHDS_Hypothesis *> &
1104 StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
1105                                          const TopoDS_Shape & aShape,
1106                                          const bool           ignoreAuxiliary)
1107 {
1108   _usedHypList.clear();
1109   _mainEdge.Nullify();
1110
1111   SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
1112   auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
1113   const bool ignoreAux = true;
1114   InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
1115
1116   // get non-auxiliary assigned to aShape
1117   int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
1118
1119   if (nbHyp == 0 && aShape.ShapeType() == TopAbs_EDGE)
1120   {
1121     // Check, if propagated from some other edge
1122     _mainEdge = StdMeshers_Propagation::GetPropagationSource( aMesh, aShape );
1123     if ( !_mainEdge.IsNull() )
1124     {
1125       // Propagation of 1D hypothesis from <aMainEdge> on this edge;
1126       // get non-auxiliary assigned to _mainEdge
1127       nbHyp = aMesh.GetHypotheses( _mainEdge, compatibleFilter, _usedHypList, true );
1128     }
1129   }
1130
1131   if (nbHyp == 0) // nothing propagated nor assigned to aShape
1132   {
1133     SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
1134     nbHyp = _usedHypList.size();
1135   }
1136   else
1137   {
1138     // get auxiliary hyps from aShape
1139     aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
1140   }
1141   if ( nbHyp > 1 && ignoreAuxiliary )
1142     _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
1143
1144   return _usedHypList;
1145 }