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