src/StdMeshers/StdMeshers_Regular_1D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : StdMeshers_Regular_1D.cxx
  25 //           Moved here from SMESH_Regular_1D.cxx
  26 //  Author : Paul RASCLE, EDF
  27 //  Module : SMESH
  28 //  $Header$
  29
  30 using namespace std;
  31
  32 #include "StdMeshers_Regular_1D.hxx"
  33 #include "StdMeshers_Distribution.hxx"
  34 #include "SMESH_Gen.hxx"
  35 #include "SMESH_Mesh.hxx"
  36 #include "SMESH_HypoFilter.hxx"
  37 #include "SMESH_subMesh.hxx"
  38
  39 #include "StdMeshers_LocalLength.hxx"
  40 #include "StdMeshers_NumberOfSegments.hxx"
  41 #include "StdMeshers_Arithmetic1D.hxx"
  42 #include "StdMeshers_StartEndLength.hxx"
  43 #include "StdMeshers_Deflection1D.hxx"
  44 #include "StdMeshers_AutomaticLength.hxx"
  45
  46 #include "SMDS_MeshElement.hxx"
  47 #include "SMDS_MeshNode.hxx"
  48 #include "SMDS_EdgePosition.hxx"
  49
  50 #include "Utils_SALOME_Exception.hxx"
  51 #include "utilities.h"
  52
  53 #include <BRep_Tool.hxx>
  54 #include <TopoDS_Edge.hxx>
  55 #include <TopoDS_Shape.hxx>
  56 #include <TopTools_ListIteratorOfListOfShape.hxx>
  57 #include <GeomAdaptor_Curve.hxx>
  58 #include <GCPnts_AbscissaPoint.hxx>
  59 #include <GCPnts_UniformAbscissa.hxx>
  60 #include <GCPnts_UniformDeflection.hxx>
  61 #include <Standard_ErrorHandler.hxx>
  62 #include <Precision.hxx>
  63 #include <Expr_GeneralExpression.hxx>
  64 #include <Expr_NamedUnknown.hxx>
  65 #include <Expr_Array1OfNamedUnknown.hxx>
  66 #include <TColStd_Array1OfReal.hxx>
  67 #include <ExprIntrp_GenExp.hxx>
  68 #include <OSD.hxx>
  69
  70 #include <string>
  71 #include <math.h>
  72
  73 using namespace std;
  74
  75 //=============================================================================
  76 /*!
  77  *
  78  */
  79 //=============================================================================
  80
  81 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
  82         SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
  83 {
  84         MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
  85         _name = "Regular_1D";
  86         _shapeType = (1 << TopAbs_EDGE);
  87
  88         _compatibleHypothesis.push_back("LocalLength");
  89         _compatibleHypothesis.push_back("NumberOfSegments");
  90         _compatibleHypothesis.push_back("StartEndLength");
  91         _compatibleHypothesis.push_back("Deflection1D");
  92         _compatibleHypothesis.push_back("Arithmetic1D");
  93         _compatibleHypothesis.push_back("AutomaticLength");
  94
  95         _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
  96 }
  97
  98 //=============================================================================
  99 /*!
 100  *
 101  */
 102 //=============================================================================
 103
 104 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
 105 {
 106 }
 107
 108 //=============================================================================
 109 /*!
 110  *
 111  */
 112 //=============================================================================
 113
 114 bool StdMeshers_Regular_1D::CheckHypothesis
 115                          (SMESH_Mesh&                          aMesh,
 116                           const TopoDS_Shape&                  aShape,
 117                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 118 {
 119   _hypType = NONE;
 120   _quadraticMesh = false;
 121
 122   const bool ignoreAuxiliaryHyps = false;
 123   const list <const SMESHDS_Hypothesis * > & hyps =
 124     GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
 125
 126   // find non-auxiliary hypothesis
 127   const SMESHDS_Hypothesis *theHyp = 0;
 128   list <const SMESHDS_Hypothesis * >::const_iterator h = hyps.begin();
 129   for ( ; h != hyps.end(); ++h ) {
 130     if ( static_cast<const SMESH_Hypothesis*>(*h)->IsAuxiliary() ) {
 131       if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 )
 132         _quadraticMesh = true;
 133     }
 134     else {
 135       if ( !theHyp )
 136         theHyp = *h; // use only the first non-auxiliary hypothesis
 137     }
 138   }
 139
 140   if ( !theHyp )
 141   {
 142     aStatus = SMESH_Hypothesis::HYP_MISSING;
 143     return false;  // can't work without a hypothesis
 144   }
 145
 146   string hypName = theHyp->GetName();
 147
 148   if (hypName == "LocalLength")
 149   {
 150     const StdMeshers_LocalLength * hyp =
 151       dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
 152     ASSERT(hyp);
 153     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength();
 154     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 155     _hypType = LOCAL_LENGTH;
 156     aStatus = SMESH_Hypothesis::HYP_OK;
 157   }
 158
 159   else if (hypName == "NumberOfSegments")
 160   {
 161     const StdMeshers_NumberOfSegments * hyp =
 162       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
 163     ASSERT(hyp);
 164     _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
 165     ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
 166     _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
 167     switch (_ivalue[ DISTR_TYPE_IND ])
 168     {
 169     case StdMeshers_NumberOfSegments::DT_Scale:
 170       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
 171       break;
 172     case StdMeshers_NumberOfSegments::DT_TabFunc:
 173       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
 174       break;
 175     case StdMeshers_NumberOfSegments::DT_ExprFunc:
 176       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
 177       break;
 178     case StdMeshers_NumberOfSegments::DT_Regular:
 179       break;
 180     default:
 181       ASSERT(0);
 182       break;
 183     }
 184     if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
 185         _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
 186         _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
 187     _hypType = NB_SEGMENTS;
 188     aStatus = SMESH_Hypothesis::HYP_OK;
 189   }
 190
 191   else if (hypName == "Arithmetic1D")
 192   {
 193     const StdMeshers_Arithmetic1D * hyp =
 194       dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
 195     ASSERT(hyp);
 196     _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
 197     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
 198     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
 199     _hypType = ARITHMETIC_1D;
 200     aStatus = SMESH_Hypothesis::HYP_OK;
 201   }
 202
 203   else if (hypName == "StartEndLength")
 204   {
 205     const StdMeshers_StartEndLength * hyp =
 206       dynamic_cast <const StdMeshers_StartEndLength * >(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 = BEG_END_LENGTH;
 212     aStatus = SMESH_Hypothesis::HYP_OK;
 213   }
 214
 215   else if (hypName == "Deflection1D")
 216   {
 217     const StdMeshers_Deflection1D * hyp =
 218       dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
 219     ASSERT(hyp);
 220     _value[ DEFLECTION_IND ] = hyp->GetDeflection();
 221     ASSERT( _value[ DEFLECTION_IND ] > 0 );
 222     _hypType = DEFLECTION;
 223     aStatus = SMESH_Hypothesis::HYP_OK;
 224   }
 225
 226   else if (hypName == "AutomaticLength")
 227   {
 228     StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
 229       (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
 230     ASSERT(hyp);
 231     _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
 232     ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
 233     _hypType = LOCAL_LENGTH;
 234     aStatus = SMESH_Hypothesis::HYP_OK;
 235   }
 236   else
 237     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 238
 239   return ( _hypType != NONE );
 240 }
 241
 242 //=======================================================================
 243 //function : compensateError
 244 //purpose  : adjust theParams so that the last segment length == an
 245 //=======================================================================
 246
 247 static void compensateError(double a1, double an,
 248                             double U1, double Un,
 249                             double             length,
 250                             GeomAdaptor_Curve& C3d,
 251                             list<double> &     theParams)
 252 {
 253   int i, nPar = theParams.size();
 254   if ( a1 + an < length && nPar > 1 )
 255   {
 256     list<double>::reverse_iterator itU = theParams.rbegin();
 257     double Ul = *itU++;
 258     // dist from the last point to the edge end <Un>, it should be equal <an>
 259     double Ln = GCPnts_AbscissaPoint::Length( C3d, Ul, Un );
 260     double dLn = an - Ln; // error of <an>
 261     if ( Abs( dLn ) <= Precision::Confusion() )
 262       return;
 263     double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
 264     double dUn = dLn * Abs( Un - U1 ) / length; // parametric error of <an>
 265     if ( dUn < 0.5 * dU ) { // last segment is a bit shorter than it should
 266       dUn = -dUn; // move the last parameter to the edge beginning
 267     }
 268     else {  // last segment is much shorter than it should -> remove the last param and
 269       theParams.pop_back(); nPar--; // move the rest points toward the edge end
 270       Ln = GCPnts_AbscissaPoint::Length( C3d, theParams.back(), Un );
 271       dUn = ( an - Ln ) * Abs( Un - U1 ) / length;
 272       if ( dUn < 0.5 * dU )
 273         dUn = -dUn;
 274     }
 275     if ( U1 > Un )
 276       dUn = -dUn;
 277     double q  = dUn / ( nPar - 1 );
 278     for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
 279       (*itU) += dUn;
 280       dUn -= q;
 281     }
 282   }
 283 }
 284
 285 static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
 286                                double length, bool theReverse,
 287                                int nbSeg, Function& func,
 288                                list<double>& theParams)
 289 {
 290   OSD::SetSignal( true );
 291
 292   if( nbSeg<=0 )
 293     return false;
 294
 295   MESSAGE( "computeParamByFunc" );
 296
 297   int nbPnt = 1 + nbSeg;
 298   vector<double> x(nbPnt, 0.);
 299
 300   if( !buildDistribution( func, 0.0, 1.0, nbSeg, x, 1E-4 ) )
 301      return false;
 302
 303   MESSAGE( "Points:\n" );
 304   char buf[1024];
 305   for( int i=0; i<=nbSeg; i++ )
 306   {
 307     sprintf(  buf, "%f\n", float(x[i] ) );
 308     MESSAGE( buf );
 309   }
 310
 311
 312
 313   // apply parameters in range [0,1] to the space of the curve
 314   double prevU = first;
 315   double sign = 1.;
 316   if (theReverse)
 317   {
 318     prevU = last;
 319     sign = -1.;
 320   }
 321   for( int i = 1; i < nbSeg; i++ )
 322   {
 323     double curvLength = length * (x[i] - x[i-1]) * sign;
 324     GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
 325     if ( !Discret.IsDone() )
 326       return false;
 327     double U = Discret.Parameter();
 328     if ( U > first && U < last )
 329       theParams.push_back( U );
 330     else
 331       return false;
 332     prevU = U;
 333   }
 334   return true;
 335 }
 336
 337 //=============================================================================
 338 /*!
 339  *
 340  */
 341 //=============================================================================
 342 bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge,
 343                                                       list<double> &     theParams,
 344                                                       const bool         theReverse) const
 345 {
 346   theParams.clear();
 347
 348   double f, l;
 349   Handle(Geom_Curve) Curve = BRep_Tool::Curve(theEdge, f, l);
 350   GeomAdaptor_Curve C3d(Curve);
 351
 352   double length = EdgeLength(theEdge);
 353
 354   switch( _hypType )
 355   {
 356   case LOCAL_LENGTH:
 357   case NB_SEGMENTS: {
 358
 359     double eltSize = 1;
 360     if ( _hypType == LOCAL_LENGTH )
 361     {
 362       // Local Length hypothesis
 363       double nbseg = ceil(length / _value[ BEG_LENGTH_IND ]); // integer sup
 364       if (nbseg <= 0)
 365         nbseg = 1;                        // degenerated edge
 366       eltSize = length / nbseg;
 367     }
 368     else
 369     {
 370       // Number Of Segments hypothesis
 371       switch (_ivalue[ DISTR_TYPE_IND ])
 372       {
 373       case StdMeshers_NumberOfSegments::DT_Scale:
 374         {
 375           int NbSegm   = _ivalue[ NB_SEGMENTS_IND ];
 376           double scale = _value[ SCALE_FACTOR_IND ];
 377
 378           if (fabs(scale - 1.0) < Precision::Confusion()) {
 379             // special case to avoid division on zero
 380             for (int i = 1; i < NbSegm; i++) {
 381               double param = f + (l - f) * i / NbSegm;
 382               theParams.push_back( param );
 383             }
 384           } else {
 385             // general case of scale distribution
 386             if ( theReverse )
 387               scale = 1.0 / scale;
 388
 389             double alpha = pow(scale, 1.0 / (NbSegm - 1));
 390             double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
 391
 392             for (int i = 1; i < NbSegm; i++) {
 393               double param = f + factor * (1.0 - pow(alpha, i));
 394               theParams.push_back( param );
 395             }
 396           }
 397           return true;
 398         }
 399         break;
 400       case StdMeshers_NumberOfSegments::DT_TabFunc:
 401         {
 402           FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
 403           return computeParamByFunc(C3d, f, l, length, theReverse,
 404                                     _ivalue[ NB_SEGMENTS_IND ], func,
 405                                     theParams);
 406         }
 407         break;
 408       case StdMeshers_NumberOfSegments::DT_ExprFunc:
 409         {
 410           FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
 411           return computeParamByFunc(C3d, f, l, length, theReverse,
 412                                     _ivalue[ NB_SEGMENTS_IND ], func,
 413                                     theParams);
 414         }
 415         break;
 416       case StdMeshers_NumberOfSegments::DT_Regular:
 417         eltSize = length / _ivalue[ NB_SEGMENTS_IND ];
 418         break;
 419       default:
 420         return false;
 421       }
 422     }
 423     GCPnts_UniformAbscissa Discret(C3d, eltSize, f, l);
 424     if ( !Discret.IsDone() )
 425       return false;
 426
 427     int NbPoints = Discret.NbPoints();
 428     for ( int i = 2; i < NbPoints; i++ )
 429     {
 430       double param = Discret.Parameter(i);
 431       theParams.push_back( param );
 432     }
 433     compensateError( eltSize, eltSize, f, l, length, C3d, theParams ); // for PAL9899
 434     return true;
 435   }
 436
 437   case BEG_END_LENGTH: {
 438
 439     // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = length
 440
 441     double a1 = _value[ BEG_LENGTH_IND ];
 442     double an = _value[ END_LENGTH_IND ];
 443     double q  = ( length - a1 ) / ( length - an );
 444
 445     double U1 = theReverse ? l : f;
 446     double Un = theReverse ? f : l;
 447     double param = U1;
 448     double eltSize = theReverse ? -a1 : a1;
 449     while ( 1 ) {
 450       // computes a point on a curve <C3d> at the distance <eltSize>
 451       // from the point of parameter <param>.
 452       GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
 453       if ( !Discret.IsDone() ) break;
 454       param = Discret.Parameter();
 455       if ( param > f && param < l )
 456         theParams.push_back( param );
 457       else
 458         break;
 459       eltSize *= q;
 460     }
 461     compensateError( a1, an, U1, Un, length, C3d, theParams );
 462     return true;
 463   }
 464
 465   case ARITHMETIC_1D: {
 466
 467     // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = length
 468
 469     double a1 = _value[ BEG_LENGTH_IND ];
 470     double an = _value[ END_LENGTH_IND ];
 471
 472     double  q = ( an - a1 ) / ( 2 *length/( a1 + an ) - 1 );
 473     int     n = int( 1 + ( an - a1 ) / q );
 474
 475     double U1 = theReverse ? l : f;
 476     double Un = theReverse ? f : l;
 477     double param = U1;
 478     double eltSize = a1;
 479     if ( theReverse ) {
 480       eltSize = -eltSize;
 481       q = -q;
 482     }
 483     while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
 484       // computes a point on a curve <C3d> at the distance <eltSize>
 485       // from the point of parameter <param>.
 486       GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
 487       if ( !Discret.IsDone() ) break;
 488       param = Discret.Parameter();
 489       if ( param > f && param < l )
 490         theParams.push_back( param );
 491       else
 492         break;
 493       eltSize += q;
 494     }
 495     compensateError( a1, an, U1, Un, length, C3d, theParams );
 496
 497     return true;
 498   }
 499
 500   case DEFLECTION: {
 501
 502     GCPnts_UniformDeflection Discret(C3d, _value[ DEFLECTION_IND ], f, l, true);
 503     if ( !Discret.IsDone() )
 504       return false;
 505
 506     int NbPoints = Discret.NbPoints();
 507     for ( int i = 2; i < NbPoints; i++ )
 508     {
 509       double param = Discret.Parameter(i);
 510       theParams.push_back( param );
 511     }
 512     return true;
 513
 514   }
 515
 516   default:;
 517   }
 518
 519   return false;
 520 }
 521
 522 //=============================================================================
 523 /*!
 524  *
 525  */
 526 //=============================================================================
 527
 528 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
 529 {
 530   MESSAGE("StdMeshers_Regular_1D::Compute");
 531
 532   if ( _hypType == NONE )
 533     return false;
 534
 535   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 536   aMesh.GetSubMesh(aShape);
 537
 538   const TopoDS_Edge & EE = TopoDS::Edge(aShape);
 539   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
 540   int shapeID = meshDS->ShapeToIndex( E );
 541
 542   double f, l;
 543   Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
 544
 545   TopoDS_Vertex VFirst, VLast;
 546   TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
 547
 548   ASSERT(!VFirst.IsNull());
 549   SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
 550   if (!lid->more())
 551   {
 552     MESSAGE (" NO NODE BUILT ON VERTEX ");
 553     return false;
 554   }
 555   const SMDS_MeshNode * idFirst = lid->next();
 556
 557   ASSERT(!VLast.IsNull());
 558   lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
 559   if (!lid->more()) {
 560     MESSAGE (" NO NODE BUILT ON VERTEX ");
 561     return false;
 562   }
 563   const SMDS_MeshNode * idLast = lid->next();
 564
 565   if (!Curve.IsNull()) {
 566     list< double > params;
 567     bool reversed = false;
 568     if ( !_mainEdge.IsNull() )
 569       reversed = aMesh.IsReversedInChain( EE, _mainEdge );
 570     try {
 571       if ( ! computeInternalParameters( E, params, reversed )) {
 572         //cout << "computeInternalParameters() failed" <<endl;
 573         return false;
 574       }
 575     }
 576     catch ( Standard_Failure ) {
 577       //cout << "computeInternalParameters() failed, Standard_Failure" <<endl;
 578       return false;
 579     }
 580
 581     // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
 582     // only internal nodes receive an edge position with param on curve
 583
 584     const SMDS_MeshNode * idPrev = idFirst;
 585     double parPrev = f;
 586     double parLast = l;
 587 //     if(reversed) {
 588 //       parPrev = l;
 589 //       parLast = f;
 590 //     }
 591
 592     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
 593       double param = *itU;
 594       gp_Pnt P = Curve->Value(param);
 595
 596       //Add the Node in the DataStructure
 597       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 598       meshDS->SetNodeOnEdge(node, shapeID, param);
 599
 600       if(_quadraticMesh) {
 601         // create medium node
 602         double prm = ( parPrev + param )/2;
 603         gp_Pnt PM = Curve->Value(prm);
 604         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 605         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 606         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 607         meshDS->SetMeshElementOnShape(edge, shapeID);
 608       }
 609       else {
 610         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 611         meshDS->SetMeshElementOnShape(edge, shapeID);
 612       }
 613
 614       idPrev = node;
 615       parPrev = param;
 616     }
 617     if(_quadraticMesh) {
 618       double prm = ( parPrev + parLast )/2;
 619       gp_Pnt PM = Curve->Value(prm);
 620       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 621       meshDS->SetNodeOnEdge(NM, shapeID, prm);
 622       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
 623       meshDS->SetMeshElementOnShape(edge, shapeID);
 624     }
 625     else {
 626       SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
 627       meshDS->SetMeshElementOnShape(edge, shapeID);
 628     }
 629   }
 630   else {
 631     // Edge is a degenerated Edge : We put n = 5 points on the edge.
 632     const int NbPoints = 5;
 633     BRep_Tool::Range(E, f, l);
 634     double du = (l - f) / (NbPoints - 1);
 635     //MESSAGE("************* Degenerated edge! *****************");
 636
 637     TopoDS_Vertex V1, V2;
 638     TopExp::Vertices(E, V1, V2);
 639     gp_Pnt P = BRep_Tool::Pnt(V1);
 640
 641     const SMDS_MeshNode * idPrev = idFirst;
 642     for (int i = 2; i < NbPoints; i++) {
 643       double param = f + (i - 1) * du;
 644       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 645       if(_quadraticMesh) {
 646         // create medium node
 647         double prm = param - du/2.;
 648         gp_Pnt PM = Curve->Value(prm);
 649         SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 650         meshDS->SetNodeOnEdge(NM, shapeID, prm);
 651         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM);
 652         meshDS->SetMeshElementOnShape(edge, shapeID);
 653       }
 654       else {
 655         SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
 656         meshDS->SetMeshElementOnShape(edge, shapeID);
 657       }
 658       meshDS->SetNodeOnEdge(node, shapeID, param);
 659       idPrev = node;
 660     }
 661     if(_quadraticMesh) {
 662       // create medium node
 663       double prm = l - du/2.;
 664       gp_Pnt PM = Curve->Value(prm);
 665       SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z());
 666       meshDS->SetNodeOnEdge(NM, shapeID, prm);
 667       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM);
 668       meshDS->SetMeshElementOnShape(edge, shapeID);
 669     }
 670     else {
 671       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
 672       meshDS->SetMeshElementOnShape(edge, shapeID);
 673     }
 674   }
 675   return true;
 676 }
 677
 678 //=============================================================================
 679 /*!
 680  *  See comments in SMESH_Algo.cxx
 681  */
 682 //=============================================================================
 683
 684 const list <const SMESHDS_Hypothesis *> &
 685 StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
 686                                          const TopoDS_Shape & aShape,
 687                                          const bool           ignoreAuxiliary)
 688 {
 689   _usedHypList.clear();
 690   _mainEdge.Nullify();
 691
 692   SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
 693   auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
 694   const bool ignoreAux = true;
 695   InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
 696
 697   // get non-auxiliary assigned to aShape
 698   int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
 699
 700   if (nbHyp == 0)
 701   {
 702     // Check, if propagated from some other edge
 703     if (aShape.ShapeType() == TopAbs_EDGE &&
 704         aMesh.IsPropagatedHypothesis(aShape, _mainEdge))
 705     {
 706       // Propagation of 1D hypothesis from <aMainEdge> on this edge;
 707       // get non-auxiliary assigned to _mainEdge
 708       nbHyp = aMesh.GetHypotheses( _mainEdge, compatibleFilter, _usedHypList, false );
 709     }
 710   }
 711
 712   if (nbHyp == 0) // nothing propagated nor assigned to aShape
 713   {
 714     SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary );
 715     nbHyp = _usedHypList.size();
 716   }
 717   else
 718   {
 719     // get auxiliary hyps from aShape
 720     aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true );
 721   }
 722   if ( nbHyp > 1 && ignoreAuxiliary )
 723     _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed
 724
 725   return _usedHypList;
 726 }
 727
 728 //=============================================================================
 729 /*!
 730  *
 731  */
 732 //=============================================================================
 733
 734 ostream & StdMeshers_Regular_1D::SaveTo(ostream & save)
 735 {
 736   return save;
 737 }
 738
 739 //=============================================================================
 740 /*!
 741  *
 742  */
 743 //=============================================================================
 744
 745 istream & StdMeshers_Regular_1D::LoadFrom(istream & load)
 746 {
 747   return load;
 748 }
 749
 750 //=============================================================================
 751 /*!
 752  *
 753  */
 754 //=============================================================================
 755
 756 ostream & operator <<(ostream & save, StdMeshers_Regular_1D & hyp)
 757 {
 758   return hyp.SaveTo( save );
 759 }
 760
 761 //=============================================================================
 762 /*!
 763  *
 764  */
 765 //=============================================================================
 766
 767 istream & operator >>(istream & load, StdMeshers_Regular_1D & hyp)
 768 {
 769   return hyp.LoadFrom( load );
 770 }