src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

   1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // File  : StdMeshers_RadialQuadrangle_1D2D.cxx
  21 // Module: SMESH
  22
  23 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
  24
  25 #include "StdMeshers_NumberOfLayers.hxx"
  26 #include "StdMeshers_LayerDistribution.hxx"
  27 #include "StdMeshers_Regular_1D.hxx"
  28 #include "StdMeshers_NumberOfSegments.hxx"
  29 #include "StdMeshers_FaceSide.hxx"
  30
  31 #include "SMDS_MeshNode.hxx"
  32 #include "SMESHDS_Mesh.hxx"
  33 #include "SMESHDS_SubMesh.hxx"
  34 #include "SMESH_Gen.hxx"
  35 #include "SMESH_HypoFilter.hxx"
  36 #include "SMESH_Mesh.hxx"
  37 #include "SMESH_MesherHelper.hxx"
  38 #include "SMESH_subMesh.hxx"
  39 #include "SMESH_subMeshEventListener.hxx"
  40 #include "SMESH_Block.hxx"
  41
  42 #include "utilities.h"
  43
  44 #include <BRepAdaptor_CompCurve.hxx>
  45 #include <BRepAdaptor_Curve.hxx>
  46 #include <BRepBuilderAPI_MakeEdge.hxx>
  47 #include <BRep_Builder.hxx>
  48 #include <BRep_Tool.hxx>
  49 #include <GCPnts_AbscissaPoint.hxx>
  50 #include <Geom2d_Line.hxx>
  51 #include <Geom2d_TrimmedCurve.hxx>
  52 #include <GeomAPI_ProjectPointOnSurf.hxx>
  53 #include <Geom_Circle.hxx>
  54 #include <Geom_Line.hxx>
  55 #include <Geom_TrimmedCurve.hxx>
  56 #include <ShapeFix_Edge.hxx>
  57 #include <TColgp_SequenceOfPnt.hxx>
  58 #include <TColgp_SequenceOfPnt2d.hxx>
  59 #include <TopExp.hxx>
  60 #include <TopExp_Explorer.hxx>
  61 #include <TopTools_ListIteratorOfListOfShape.hxx>
  62 #include <TopoDS.hxx>
  63
  64
  65 using namespace std;
  66
  67 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  68 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  69
  70
  71 //=======================================================================
  72 //function : StdMeshers_RadialQuadrangle_1D2D
  73 //purpose  :
  74 //=======================================================================
  75
  76 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int        hypId,
  77                                                                    SMESH_Gen* gen)
  78   :StdMeshers_Quadrangle_2D( hypId, gen )
  79 {
  80   _name = "RadialQuadrangle_1D2D";
  81   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
  82
  83   _compatibleHypothesis.push_back("LayerDistribution2D");
  84   _compatibleHypothesis.push_back("NumberOfLayers2D");
  85   _requireDiscreteBoundary = false;
  86   _supportSubmeshes = true;
  87   _neededLowerHyps[ 1 ] = true;  // suppress warning on hiding a global 1D algo
  88
  89   myNbLayerHypo      = 0;
  90   myDistributionHypo = 0;
  91 }
  92
  93
  94 //================================================================================
  95 /*!
  96  * \brief Destructor
  97  */
  98 //================================================================================
  99
 100 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
 101 {}
 102
 103
 104 //=======================================================================
 105 //function : CheckHypothesis
 106 //purpose  :
 107 //=======================================================================
 108
 109 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
 110                            (SMESH_Mesh&                          aMesh,
 111                             const TopoDS_Shape&                  aShape,
 112                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
 113 {
 114   // check aShape
 115   myNbLayerHypo = 0;
 116   myDistributionHypo = 0;
 117
 118   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 119
 120   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 121   if ( hyps.size() == 0 ) {
 122     aStatus = SMESH_Hypothesis::HYP_OK;
 123     return true;  // can work with no hypothesis
 124   }
 125
 126   if ( hyps.size() > 1 ) {
 127     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 128     return false;
 129   }
 130
 131   const SMESHDS_Hypothesis *theHyp = hyps.front();
 132
 133   string hypName = theHyp->GetName();
 134
 135   if (hypName == "NumberOfLayers2D") {
 136     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 137     aStatus = SMESH_Hypothesis::HYP_OK;
 138     return true;
 139   }
 140   if (hypName == "LayerDistribution2D") {
 141     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 142     aStatus = SMESH_Hypothesis::HYP_OK;
 143     return true;
 144   }
 145   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 146   return true;
 147 }
 148
 149 namespace
 150 {
 151   // ------------------------------------------------------------------------------
 152   /*!
 153    * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
 154    */
 155   class TEdgeMarker : public SMESH_subMeshEventListener
 156   {
 157     TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
 158                                               "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
 159   public:
 160     //!<  Return static listener
 161     static SMESH_subMeshEventListener* getListener()
 162     {
 163       static TEdgeMarker theEdgeMarker;
 164       return &theEdgeMarker;
 165     }
 166     //! Clear edge sumbesh if something happens on face
 167     void ProcessEvent(const int          /*event*/,
 168                       const int          eventType,
 169                       SMESH_subMesh*     /*faceSubMesh*/,
 170                       EventListenerData* edgesHolder,
 171                       const SMESH_Hypothesis*  /*hyp*/)
 172     {
 173       if ( edgesHolder && eventType == SMESH_subMesh::ALGO_EVENT)
 174       {
 175         std::list<SMESH_subMesh*>::iterator smIt = edgesHolder->mySubMeshes.begin();
 176         for ( ; smIt != edgesHolder->mySubMeshes.end(); ++smIt )
 177         {
 178           SMESH_subMesh* edgeSM = *smIt;
 179           edgeSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
 180         }
 181       }
 182     }
 183     //! Store edge SMESH_subMesh'es computed by the algo
 184     static void markEdge( const TopoDS_Edge& edge, SMESH_subMesh* faceSM )
 185     {
 186       if ( SMESH_subMesh* edgeSM = faceSM->GetFather()->GetSubMeshContaining( edge ))
 187       {
 188         EventListenerData* edgesHolder = faceSM->GetEventListenerData( getListener() );
 189         if ( edgesHolder )
 190         {
 191           std::list<SMESH_subMesh*>::iterator smIt = std::find( edgesHolder->mySubMeshes.begin(),
 192                                                                 edgesHolder->mySubMeshes.end(),
 193                                                                 edgeSM );
 194           if ( smIt == edgesHolder->mySubMeshes.end() )
 195             edgesHolder->mySubMeshes.push_back( edgeSM );
 196         }
 197         else
 198         {
 199           edgesHolder = SMESH_subMeshEventListenerData::MakeData( edgeSM );
 200           faceSM->SetEventListener( TEdgeMarker::getListener(), edgesHolder, faceSM );
 201         }
 202       }
 203     }
 204   };
 205
 206   //================================================================================
 207   /*!
 208    * \brief Return sides of the face connected in the order: aCircEdge, aLinEdge1, aLinEdge2
 209    *  \retval int - nb of sides
 210    */
 211   //================================================================================
 212
 213   int analyseFace(const TopoDS_Shape&     aShape,
 214                   SMESH_Mesh*             aMesh,
 215                   StdMeshers_FaceSidePtr& aCircSide,
 216                   StdMeshers_FaceSidePtr& aLinSide1,
 217                   StdMeshers_FaceSidePtr& aLinSide2,
 218                   SMESH_MesherHelper*     helper)
 219   {
 220     const TopoDS_Face& face = TopoDS::Face( aShape );
 221     aCircSide.reset(); aLinSide1.reset(); aLinSide2.reset();
 222
 223     list< TopoDS_Edge > edges;
 224     list< int > nbEdgesInWire;
 225     int nbWire = SMESH_Block::GetOrderedEdges ( face, edges, nbEdgesInWire );
 226     if ( nbWire > 2 || nbEdgesInWire.front() < 1 ) return 0;
 227
 228     // remove degenerated EDGEs
 229     TopTools_MapOfShape degenVV;
 230     list<TopoDS_Edge>::iterator edge = edges.begin();
 231     while ( edge != edges.end() )
 232       if ( SMESH_Algo::isDegenerated( *edge ))
 233       {
 234         degenVV.Add( SMESH_MesherHelper::IthVertex( 0, *edge ));
 235         degenVV.Add( SMESH_MesherHelper::IthVertex( 1, *edge ));
 236         edge = edges.erase( edge );
 237       }
 238       else
 239       {
 240         ++edge;
 241       }
 242     int nbEdges = edges.size();
 243
 244     // find VERTEXes between continues EDGEs
 245     TopTools_MapOfShape contVV;
 246     if ( nbEdges > 1 )
 247     {
 248       TopoDS_Edge ePrev = edges.back();
 249       for ( edge = edges.begin(); edge != edges.end(); ++edge )
 250       {
 251         if ( SMESH_Algo::IsContinuous( ePrev, *edge ))
 252           contVV.Add( SMESH_MesherHelper::IthVertex( 0, *edge ));
 253         ePrev = *edge;
 254       }
 255     }
 256     // make edges start from a non-continues VERTEX
 257     if ( 1 < contVV.Extent() && contVV.Extent() < nbEdges )
 258     {
 259       while ( contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
 260         edges.splice( edges.end(), edges, edges.begin() );
 261     }
 262
 263     // make face sides
 264     TSideVector sides;
 265     while ( !edges.empty() )
 266     {
 267       list< TopoDS_Edge > sideEdges;
 268       sideEdges.splice( sideEdges.end(), edges, edges.begin() );
 269       while ( !edges.empty() &&
 270               contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
 271         sideEdges.splice( sideEdges.end(), edges, edges.begin() );
 272
 273       StdMeshers_FaceSidePtr side;
 274       if ( aMesh )
 275         side = StdMeshers_FaceSide::New( face, sideEdges, aMesh,
 276                                          /*isFwd=*/true, /*skipMedium=*/ true, helper );
 277       sides.push_back( side );
 278     }
 279
 280     if ( !aMesh ) // call from IsApplicable()
 281       return sides.size();
 282
 283     if ( sides.size() > 3 )
 284       return sides.size();
 285
 286     if ( nbWire == 2 && (( sides.size() != 2 ) ||
 287                          ( sides[0]->IsClosed() && sides[1]->IsClosed() ) ||
 288                          ( !sides[0]->IsClosed() && !sides[1]->IsClosed() )))
 289       return -1;
 290
 291     // detect an elliptic side
 292
 293     if ( sides.size() == 1 )
 294     {
 295       aCircSide = sides[0];
 296       return sides.size();
 297     }
 298
 299     // sort sides by deviation from a straight line
 300     multimap< double, int > deviation2sideInd;
 301     const double nbSamples = 7;
 302     for ( size_t iS = 0; iS < sides.size(); ++iS )
 303     {
 304       gp_Pnt pf = BRep_Tool::Pnt( sides[iS]->FirstVertex() );
 305       gp_Pnt pl = BRep_Tool::Pnt( sides[iS]->LastVertex() );
 306       gp_Vec v1( pf, pl );
 307       double v1Len = v1.Magnitude();
 308       if ( v1Len < std::numeric_limits< double >::min() )
 309       {
 310         deviation2sideInd.insert( make_pair( sides[iS]->Length(), iS )); // the side seems closed
 311         continue;
 312       }
 313       double devia = 0;
 314       for ( int i = 0; i < nbSamples; ++i )
 315       {
 316         gp_Pnt pi( sides[iS]->Value3d(( i + 1 ) / nbSamples ));
 317         gp_Vec vi( pf, pi );
 318         double h = 0.5 * v1.Crossed( vi ).Magnitude() / v1Len;
 319         devia = Max( devia, h );
 320       }
 321       deviation2sideInd.insert( make_pair( devia, iS ));
 322     }
 323     double maxDevi = deviation2sideInd.rbegin()->first;
 324     if ( maxDevi < 1e-7 && sides.size() == 3 )
 325     {
 326       // a triangle FACE; use a side with the most outstanding length as an elliptic one
 327       deviation2sideInd.clear();
 328       multimap< double, int > len2sideInd;
 329       for ( size_t iS = 0; iS < sides.size(); ++iS )
 330         len2sideInd.insert( make_pair( sides[iS]->Length(), iS ));
 331
 332       multimap< double, int >::iterator l2i = len2sideInd.begin();
 333       double len0 = l2i->first;
 334       double len1 = (++l2i)->first;
 335       double len2 = (++l2i)->first;
 336       if ( len1 - len0 > len2 - len1 )
 337         deviation2sideInd.insert( std::make_pair( 0., len2sideInd.begin()->second ));
 338       else
 339         deviation2sideInd.insert( std::make_pair( 0., len2sideInd.rbegin()->second ));
 340     }
 341
 342     double minDevi = deviation2sideInd.begin()->first;
 343     int   iMinCurv = deviation2sideInd.begin()->second;
 344     if ( sides.size() == 3 && degenVV.Size() == 1 &&
 345          minDevi / sides[ iMinCurv ]->Length() > 1e-3 )
 346     {
 347       // a triangle with curved sides and a degenerated EDGE (IPAL54585);
 348       // use a side opposite to the degenerated EDGE as an elliptic one
 349       for ( size_t iS = 0; iS < sides.size(); ++iS )
 350         if ( degenVV.Contains( sides[ iS ]->FirstVertex() ))
 351         {
 352           deviation2sideInd.clear();
 353           deviation2sideInd.insert( std::make_pair( 0.,( iS + 1 ) % sides.size() ));
 354           break;
 355         }
 356     }
 357
 358     int iCirc = deviation2sideInd.rbegin()->second;
 359     aCircSide = sides[ iCirc ];
 360     aLinSide1 = sides[( iCirc + 1 ) % sides.size() ];
 361     if ( sides.size() > 2 )
 362     {
 363       aLinSide2 = sides[( iCirc + 2 ) % sides.size() ];
 364       aLinSide2->Reverse(); // to be "parallel" to aLinSide1
 365     }
 366
 367     if (( nbWire == 2 && aLinSide1 ) &&
 368         ( aLinSide1->Edge(0).Orientation() == TopAbs_INTERNAL ) &&
 369         ( aCircSide->IsClosed() ))
 370     {
 371       // assure that aCircSide starts at aLinSide1
 372       TopoDS_Vertex v0 = aLinSide1->FirstVertex();
 373       TopoDS_Vertex v1 = aLinSide1->LastVertex();
 374       if ( ! aCircSide->FirstVertex().IsSame( v0 ) &&
 375            ! aCircSide->FirstVertex().IsSame( v1 ))
 376       {
 377         int iE = 0;
 378         for ( ; iE < aCircSide->NbEdges(); ++iE )
 379           if ( aCircSide->FirstVertex(iE).IsSame( v0 ) ||
 380                aCircSide->FirstVertex(iE).IsSame( v1 ))
 381             break;
 382         if ( iE == aCircSide->NbEdges() )
 383           return -2;
 384
 385         edges.clear();
 386         for ( int i = 0; i < aCircSide->NbEdges(); ++i, ++iE )
 387           edges.push_back( aCircSide->Edge( iE % aCircSide->NbEdges() ));
 388
 389         aCircSide = StdMeshers_FaceSide::New( face, edges, aMesh,
 390                                               /*isFwd=*/true, /*skipMedium=*/ true, helper );
 391       }
 392     }
 393
 394     return sides.size();
 395   }
 396
 397   //================================================================================
 398   /*!
 399    * \brief Checks if the common vertex between LinSide's lies inside the circle
 400    *  and not outside
 401    *  \return bool - false if there are 3 EDGEs and the corner is outside
 402    */
 403   //================================================================================
 404
 405   bool isCornerInsideCircle(const StdMeshers_FaceSidePtr& /*CircSide*/,
 406                             const StdMeshers_FaceSidePtr& /*LinSide1*/,
 407                             const StdMeshers_FaceSidePtr& /*LinSide2*/)
 408   {
 409     // if ( CircSide && LinSide1 && LinSide2 )
 410     // {
 411     //   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircSide ));
 412     //   TopoDS_Vertex aCommonV;
 413     //   if ( !aCirc.IsNull() &&
 414     //        TopExp::CommonVertex( LinSide1, LinSide2, aCommonV ))
 415     //   {
 416     //     gp_Pnt aCommonP = BRep_Tool::Pnt( aCommonV );
 417     //     gp_Pnt  aCenter = aCirc->Location();
 418     //     double     dist = aCenter.Distance( aCommonP );
 419     //     return dist < 0.1 * aCirc->Radius();
 420     //   }
 421     // }
 422     return true;
 423   }
 424
 425   //================================================================================
 426   /*!
 427    * \brief Create an EDGE connecting the ellipse center with the most distant point
 428    *        of the ellipse.
 429    *  \param [in] circSide - the elliptic side
 430   *  \param [in] face - the FACE
 431   *  \param [out] circNode - a node on circSide most distant from the center
 432   *  \return TopoDS_Edge - the create EDGE
 433   */
 434   //================================================================================
 435
 436   TopoDS_Edge makeEdgeToCenter( StdMeshers_FaceSidePtr& circSide,
 437                                 const TopoDS_Face&      face,
 438                                 const SMDS_MeshNode*&   circNode)
 439   {
 440     // find the center and a point most distant from it
 441
 442     double maxDist = 0, normPar = 0;
 443     gp_XY uv1, uv2;
 444     for ( int i = 0; i < 32; ++i )
 445     {
 446       double    u = 0.5 * i / 32.;
 447       gp_Pnt2d p1 = circSide->Value2d( u );
 448       gp_Pnt2d p2 = circSide->Value2d( u + 0.5 );
 449       double dist = p1.SquareDistance( p2 );
 450       if ( dist > maxDist )
 451       {
 452         maxDist = dist;
 453         uv1 = p1.XY();
 454         uv2 = p2.XY();
 455         normPar = u;
 456       }
 457     }
 458     gp_XY center = 0.5 * ( uv1 + uv2 );
 459     double   len = 0.5 * Sqrt( maxDist );
 460     bool  isCirc = ( Abs( len - circSide->Value2d( 0 ).Distance( center )) < 1e-3 * len );
 461
 462     // find a node closest to the most distant point
 463
 464     size_t iDist = 0;
 465     const UVPtStructVec& circNodes = circSide->GetUVPtStruct();
 466     if ( !isCirc )
 467     {
 468       double minDist = 1e100;
 469       for ( size_t i = 0; i <= circNodes.size(); ++i )
 470       {
 471         double dist = Abs( circNodes[i].normParam - normPar );
 472         if ( dist < minDist )
 473         {
 474           iDist = i;
 475           minDist = dist;
 476         }
 477       }
 478     }
 479     circNode = circNodes[iDist].node;
 480     uv1      = circNodes[iDist].UV();
 481     len      = ( uv1 - center ).Modulus();
 482
 483     // make the EDGE between the most distant point and the center
 484
 485     Handle(Geom2d_Line) line = new Geom2d_Line( uv1, gp_Dir2d( center - uv1 ) );
 486     Handle(Geom2d_Curve) pcu = new Geom2d_TrimmedCurve( line, 0, len );
 487
 488     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
 489     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pcu, surface, 0, len );
 490
 491     BRep_Builder().UpdateEdge( edge, pcu, face, 1e-7 );
 492     ShapeFix_Edge().FixAddCurve3d( edge );
 493
 494
 495     // assure that circSide starts at circNode
 496     if ( iDist != 0  &&  iDist != circNodes.size()-1 )
 497     {
 498       // create a new circSide
 499       UVPtStructVec nodesNew;
 500       nodesNew.reserve( circNodes.size() );
 501       nodesNew.insert( nodesNew.end(), circNodes.begin() + iDist, circNodes.end() );
 502       nodesNew.insert( nodesNew.end(), circNodes.begin() + 1, circNodes.begin() + iDist + 1 );
 503       circSide = StdMeshers_FaceSide::New( nodesNew );
 504     }
 505
 506     return edge;
 507   }
 508
 509   //================================================================================
 510   /*!
 511    * \brief Set nodes existing on a linSide to UVPtStructVec and create missing nodes
 512    *        corresponding to layerPositions
 513    */
 514   //================================================================================
 515
 516   void makeMissingMesh( StdMeshers_FaceSidePtr& linSide,
 517                         UVPtStructVec&          nodes,
 518                         const vector< double >& layerPositions,
 519                         SMESH_MesherHelper*     helper )
 520   {
 521     // tolerance to compare normParam
 522     double tol = 1e100;
 523     for ( size_t i = 1; i < layerPositions.size(); ++i )
 524       tol = Min( tol, layerPositions[i] - layerPositions[i-1] );
 525     tol *= 0.05;
 526
 527     // merge existing nodes with layerPositions into UVPtStructVec
 528     // ------------------------------------------------------------
 529
 530     const UVPtStructVec& exiNodes = linSide->GetUVPtStruct();
 531     nodes.clear();
 532     nodes.reserve( layerPositions.size() + exiNodes.size() );
 533     vector< double >::const_iterator pos = layerPositions.begin(), posEnd = layerPositions.end();
 534     for ( size_t i = 0; i < exiNodes.size(); ++i )
 535     {
 536       switch ( exiNodes[i].node->GetPosition()->GetTypeOfPosition() )
 537       {
 538       case SMDS_TOP_VERTEX:
 539       {
 540         // allocate UVPtStruct's for non-existing nodes
 541         while ( pos != posEnd && *pos < exiNodes[i].normParam - tol )
 542         {
 543           UVPtStruct uvPS;
 544           uvPS.normParam = *pos++;
 545           nodes.push_back( uvPS );
 546         }
 547         // save existing node on a VERTEX
 548         nodes.push_back( exiNodes[i] );
 549         break;
 550       }
 551       case SMDS_TOP_EDGE:
 552       {
 553         // save existing nodes on an EDGE
 554         while ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 1 )
 555         {
 556           nodes.push_back( exiNodes[i++] );
 557         }
 558         // save existing node on a VERTEX
 559         if ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 0 )
 560         {
 561           nodes.push_back( exiNodes[i] );
 562         }
 563         break;
 564       }
 565       default:;
 566       }
 567
 568       // skip layer positions covered by saved nodes
 569       while ( pos != posEnd && *pos < nodes.back().normParam + tol )
 570       {
 571         ++pos;
 572       }
 573     }
 574     // allocate UVPtStruct's for the rest non-existing nodes
 575     while ( pos != posEnd )
 576     {
 577       UVPtStruct uvPS;
 578       uvPS.normParam = *pos++;
 579       nodes.push_back( uvPS );
 580     }
 581
 582     // create missing nodes
 583     // ---------------------
 584
 585     SMESHDS_Mesh * meshDS = helper->GetMeshDS();
 586     const TopoDS_Face& F  = TopoDS::Face( helper->GetSubShape() );
 587     Handle(ShapeAnalysis_Surface) surface = helper->GetSurface( F );
 588
 589     helper->SetElementsOnShape( false ); // we create nodes on EDGEs, not on the FACE
 590
 591     for ( size_t i = 0; i < nodes.size(); ++i )
 592     {
 593       if ( nodes[ i ].node ) continue;
 594
 595       gp_Pnt2d uv = linSide->Value2d( nodes[i].normParam );
 596       gp_Pnt  xyz = surface->Value( uv.X(), uv.Y() );
 597
 598       nodes[ i ].SetUV( uv.XY() );
 599       nodes[ i ].node = helper->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
 600     }
 601
 602     // set nodes on VERTEXes
 603     for ( int iE = 0; iE < linSide->NbEdges(); ++iE )
 604     {
 605       TopoDS_Vertex v = linSide->LastVertex( iE );
 606       if ( SMESH_Algo::VertexNode( v, meshDS ))
 607         continue;
 608
 609       double normPar = linSide->LastParameter( iE );
 610       size_t i = 0;
 611       while ( nodes[ i ].normParam < normPar )
 612         ++i;
 613       if (( nodes[ i ].normParam - normPar ) > ( normPar - nodes[ i-1 ].normParam ))
 614         --i;
 615       meshDS->SetNodeOnVertex( nodes[ i ].node, v );
 616     }
 617
 618     // set nodes on EDGEs
 619     int edgeID;
 620     for ( size_t i = 0; i < nodes.size(); ++i )
 621     {
 622       if ( nodes[ i ].node->getshapeId() > 0 ) continue;
 623
 624       double u = linSide->Parameter( nodes[i].normParam, edgeID );
 625       meshDS->SetNodeOnEdge( nodes[ i ].node, edgeID, u );
 626     }
 627
 628     // create segments
 629     for ( size_t i = 1; i < nodes.size(); ++i )
 630     {
 631       if ( meshDS->FindEdge( nodes[i].node, nodes[i-1].node )) continue;
 632
 633       const SMDS_MeshElement* seg = helper->AddEdge( nodes[i].node, nodes[i-1].node );
 634
 635       double normParam = 0.5 * ( nodes[i].normParam + nodes[i-1].normParam );
 636       edgeID = linSide->EdgeID( linSide->EdgeIndex( normParam ));
 637       meshDS->SetMeshElementOnShape( seg, edgeID );
 638     }
 639
 640     helper->SetElementsOnShape( true );
 641
 642   } // end makeMissingMesh()
 643
 644 //================================================================================
 645 //================================================================================
 646 /*!
 647  * \brief Class computing layers distribution using data of
 648  *        StdMeshers_LayerDistribution hypothesis
 649  */
 650 //================================================================================
 651 //================================================================================
 652
 653 class TNodeDistributor: public StdMeshers_Regular_1D
 654 {
 655   list <const SMESHDS_Hypothesis *> myUsedHyps;
 656 public:
 657   // -----------------------------------------------------------------------------
 658   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 659   {
 660     const int myID = -1001;
 661     TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
 662     if ( !myHyp )
 663       myHyp = new TNodeDistributor( myID, aMesh.GetGen() );
 664     return myHyp;
 665   }
 666   // -----------------------------------------------------------------------------
 667   //! Computes distribution of nodes on a straight line ending at pIn and pOut
 668   bool Compute( vector< double > &      positions,
 669                 const TopoDS_Edge&      edge,
 670                 Adaptor3d_Curve&        curve,
 671                 double                  f,
 672                 double                  l,
 673                 SMESH_Mesh&             mesh,
 674                 const SMESH_Hypothesis* hyp1d)
 675   {
 676     if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
 677
 678     myUsedHyps.clear();
 679     myUsedHyps.push_back( hyp1d );
 680
 681     SMESH_Hypothesis::Hypothesis_Status aStatus;
 682     if ( !StdMeshers_Regular_1D::CheckHypothesis( mesh, edge, aStatus ))
 683       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 684                     "with LayerDistribution hypothesis");
 685
 686     double len = GCPnts_AbscissaPoint::Length( curve, f, l );
 687
 688     list< double > params;
 689     if ( !StdMeshers_Regular_1D::computeInternalParameters( mesh, curve, len, f, l, params, false ))
 690       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 691
 692     params.push_front( f );
 693     params.push_back ( l );
 694     positions.clear();
 695     positions.reserve( params.size() );
 696     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 697       positions.push_back(( *itU - f ) / ( l - f ));
 698     return true;
 699   }
 700   // -----------------------------------------------------------------------------
 701   //! Make mesh on an edge using assigned 1d hyp or default nb of segments
 702   bool ComputeCircularEdge( SMESH_Mesh&                   aMesh,
 703                             const StdMeshers_FaceSidePtr& aSide )
 704   {
 705     bool ok = true;
 706     for ( int i = 0; i < aSide->NbEdges(); ++i )
 707     {
 708       const TopoDS_Edge& edge = aSide->Edge( i );
 709       _gen->Compute( aMesh, edge );
 710       SMESH_subMesh *sm = aMesh.GetSubMesh( edge );
 711       if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
 712       {
 713         // find any 1d hyp assigned (there can be a hyp w/o algo)
 714         myUsedHyps = SMESH_Algo::GetUsedHypothesis( aMesh, edge, /*ignoreAux=*/true );
 715         Hypothesis_Status aStatus;
 716         if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 717         {
 718           // no valid 1d hyp assigned, use default nb of segments
 719           _hypType                    = NB_SEGMENTS;
 720           _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 721           _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 722         }
 723         ok &= StdMeshers_Regular_1D::Compute( aMesh, edge );
 724       }
 725     }
 726     return ok;
 727   }
 728   // -----------------------------------------------------------------------------
 729   //! Make mesh on an edge using assigned 1d hyp or default nb of segments
 730   bool EvaluateCircularEdge(SMESH_Mesh&                  aMesh,
 731                             const StdMeshers_FaceSidePtr aSide,
 732                             MapShapeNbElems&             aResMap)
 733   {
 734     bool ok = true;
 735     for ( int i = 0; i < aSide->NbEdges(); ++i )
 736     {
 737       const TopoDS_Edge& anEdge = aSide->Edge( i );
 738       _gen->Evaluate( aMesh, anEdge, aResMap );
 739       if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
 740         continue;
 741
 742       // find any 1d hyp assigned
 743       myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
 744       Hypothesis_Status aStatus;
 745       if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
 746       {
 747         // no valid 1d hyp assigned, use default nb of segments
 748         _hypType                    = NB_SEGMENTS;
 749         _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 750         _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 751       }
 752       ok &= StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
 753     }
 754     return ok;
 755   }
 756 protected:
 757   // -----------------------------------------------------------------------------
 758   TNodeDistributor( int hypId, SMESH_Gen* gen)
 759     : StdMeshers_Regular_1D( hypId, gen)
 760   {
 761   }
 762   // -----------------------------------------------------------------------------
 763   virtual const list <const SMESHDS_Hypothesis *> &
 764     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 765   {
 766     return myUsedHyps;
 767   }
 768   // -----------------------------------------------------------------------------
 769 };
 770 }
 771
 772 //=======================================================================
 773 /*!
 774  * \brief Allow algo to do something after persistent restoration
 775  * \param subMesh - restored submesh
 776  *
 777  * call TEdgeMarker::markEdge()
 778  */
 779 //=======================================================================
 780
 781 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
 782 {
 783   if ( !faceSubMesh->IsEmpty() )
 784     SetEventListener( faceSubMesh );
 785 }
 786
 787 //=======================================================================
 788 /*!
 789  * \brief Sets event listener to a submesh
 790  * \param subMesh - submesh where algo is set
 791  *
 792  * This method is called when a submesh gets HYP_OK algo_state.
 793  */
 794 //=======================================================================
 795
 796 void StdMeshers_RadialQuadrangle_1D2D::SetEventListener(SMESH_subMesh* faceSubMesh)
 797 {
 798   for ( TopExp_Explorer e( faceSubMesh->GetSubShape(), TopAbs_EDGE ); e.More(); e.Next() )
 799   {
 800     TEdgeMarker::markEdge( TopoDS::Edge( e.Current() ), faceSubMesh );
 801   }
 802 }
 803
 804 //=======================================================================
 805 //function : Compute
 806 //purpose  :
 807 //=======================================================================
 808
 809 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 810                                                const TopoDS_Shape& aShape)
 811 {
 812   SMESH_MesherHelper helper( aMesh );
 813   StdMeshers_Quadrangle_2D::myHelper     = & helper;
 814   StdMeshers_Quadrangle_2D::myNeedSmooth = false;
 815   StdMeshers_Quadrangle_2D::myCheckOri   = false;
 816   StdMeshers_Quadrangle_2D::myQuadList.clear();
 817
 818   myHelper->SetSubShape( aShape );
 819   myHelper->SetElementsOnShape( true );
 820
 821   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
 822   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
 823   if( nbSides > 3 || nbSides < 1 )
 824     return error("The face must be a full ellipse or a part of ellipse (i.e. the number "
 825                  "of edges is less or equal to 3 and one of them is an ellipse curve)");
 826
 827   // get not yet computed EDGEs
 828   list< TopoDS_Edge > emptyEdges;
 829   for ( TopExp_Explorer e( aShape, TopAbs_EDGE ); e.More(); e.Next() )
 830   {
 831     if ( aMesh.GetSubMesh( e.Current() )->IsEmpty() )
 832       emptyEdges.push_back( TopoDS::Edge( e.Current() ));
 833   }
 834
 835   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
 836
 837   if ( !algo1d->ComputeCircularEdge( aMesh, circSide ))
 838     return error( algo1d->GetComputeError() );
 839
 840
 841   TopoDS_Face F = TopoDS::Face(aShape);
 842   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 843
 844   myHelper->IsQuadraticSubMesh( aShape );
 845
 846   vector< double > layerPositions; // [0,1]
 847
 848   const SMDS_MeshNode* centerNode = 0;
 849   gp_Pnt2d             centerUV(0,0);
 850
 851   // ------------------------------------------------------------------------------------------
 852   if ( nbSides == 1 ) // full ellipse
 853   {
 854     const SMDS_MeshNode* circNode;
 855     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
 856
 857     StdMeshers_FaceSidePtr tmpSide =
 858       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
 859
 860     if ( !computeLayerPositions( tmpSide, layerPositions ))
 861       return false;
 862
 863     UVPtStructVec nodes( layerPositions.size() );
 864     nodes[0].node = circNode;
 865     for ( size_t i = 0; i < layerPositions.size(); ++i )
 866     {
 867       gp_Pnt2d uv = tmpSide->Value2d( layerPositions[i] );
 868       gp_Pnt  xyz = S->Value( uv.X(), uv.Y() );
 869
 870       nodes[ i ].SetUV( uv.XY() );
 871       nodes[ i ].normParam = layerPositions[i];
 872       if ( i )
 873         nodes[ i ].node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, uv.X(), uv.Y() );
 874     }
 875
 876     linSide1 = StdMeshers_FaceSide::New( nodes );
 877     linSide2 = StdMeshers_FaceSide::New( nodes );
 878
 879     centerNode = nodes.back().node;
 880     centerUV   = nodes.back().UV();
 881   }
 882   // ------------------------------------------------------------------------------------------
 883   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
 884   {
 885     // full ellipse with an internal radial side
 886
 887     // eliminate INTERNAL orientation
 888     list< TopoDS_Edge > edges;
 889     for ( int iE = 0; iE < linSide1->NbEdges(); ++iE )
 890     {
 891       edges.push_back( linSide1->Edge( iE ));
 892       edges.back().Orientation( TopAbs_FORWARD );
 893     }
 894
 895     // orient the internal side
 896     bool isVIn0Shared = false;
 897     TopoDS_Vertex vIn0 = myHelper->IthVertex( 0, edges.front() );
 898     for ( int iE = 0; iE < circSide->NbEdges() && !isVIn0Shared; ++iE )
 899       isVIn0Shared = vIn0.IsSame( circSide->FirstVertex( iE ));
 900
 901     linSide1 = StdMeshers_FaceSide::New( F, edges, &aMesh,
 902                                          /*isFrw=*/isVIn0Shared, /*skipMedium=*/true, myHelper );
 903
 904     int nbMeshedEdges;
 905     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges ))
 906       return false;
 907
 908     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 909     UVPtStructVec nodes;
 910     if ( nbMeshedEdges != linSide1->NbEdges() )
 911       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 912     else
 913       nodes = linSide1->GetUVPtStruct();
 914
 915     linSide1 = StdMeshers_FaceSide::New( nodes );
 916     linSide2 = StdMeshers_FaceSide::New( nodes );
 917
 918     centerNode = nodes.back().node;
 919     centerUV   = nodes.back().UV();
 920   }
 921   // ------------------------------------------------------------------------------------------
 922   else if ( nbSides == 2 )
 923   {
 924     // find positions of layers for the first half of linSide1
 925     int nbMeshedEdges;
 926     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges, /*useHalf=*/true ))
 927       return false;
 928
 929     // make positions for the whole linSide1
 930     for ( size_t i = 0; i < layerPositions.size(); ++i )
 931     {
 932       layerPositions[i] *= 0.5;
 933     }
 934     layerPositions.reserve( layerPositions.size() * 2 );
 935     for ( int nb = layerPositions.size()-1; nb > 0; --nb )
 936       layerPositions.push_back( layerPositions.back() + layerPositions[nb] - layerPositions[nb-1] );
 937
 938     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 939     UVPtStructVec nodes;
 940     if ( nbMeshedEdges != linSide1->NbEdges() )
 941       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 942     else
 943       nodes = linSide1->GetUVPtStruct();
 944
 945     // find a central node
 946     size_t i = 0;
 947     while ( nodes[ i ].normParam < 0.5 ) ++i;
 948     if (( nodes[ i ].normParam - 0.5 ) > ( 0.5 - nodes[ i-1 ].normParam )) --i;
 949
 950     // distribute nodes between two linear sides
 951     UVPtStructVec nodes2( nodes.rbegin(), nodes.rbegin() + nodes.size() - i );
 952     nodes.resize( i + 1 );
 953
 954     linSide1 = StdMeshers_FaceSide::New( nodes );
 955     linSide2 = StdMeshers_FaceSide::New( nodes2 );
 956
 957     centerNode = nodes.back().node;
 958     centerUV   = nodes.back().UV();
 959   }
 960   // ------------------------------------------------------------------------------------------
 961   else // nbSides == 3
 962   {
 963     // one curve must be a part of ellipse and 2 other curves must be segments of line
 964
 965     int nbMeshedEdges1, nbMeshedEdges2;
 966     vector< double > layerPositions2;
 967     bool ok1 = computeLayerPositions( linSide1, layerPositions,  &nbMeshedEdges1 );
 968     bool ok2 = computeLayerPositions( linSide2, layerPositions2, &nbMeshedEdges2 );
 969     if ( !ok1 && !ok2 )
 970       return false;
 971
 972     bool linSide1Computed = ( nbMeshedEdges1 == linSide1->NbEdges() );
 973     bool linSide2Computed = ( nbMeshedEdges2 == linSide2->NbEdges() );
 974
 975     UVPtStructVec nodes;
 976
 977     if ( linSide1Computed && !linSide2Computed )
 978     {
 979       // use layer positions of linSide1 to mesh linSide2
 980       makeMissingMesh( linSide2, nodes, layerPositions, myHelper );
 981       linSide2 = StdMeshers_FaceSide::New( nodes );
 982     }
 983     else if ( linSide2Computed && !linSide1Computed )
 984     {
 985       // use layer positions of linSide2 to mesh linSide1
 986       makeMissingMesh( linSide1, nodes, layerPositions2, myHelper );
 987       linSide1 = StdMeshers_FaceSide::New( nodes );
 988     }
 989     else if ( !linSide2Computed && !linSide1Computed )
 990     {
 991       // use layer positions of a longer side to mesh the shorter side
 992       vector< double >& lp =
 993         ( linSide1->Length() > linSide2->Length() ) ? layerPositions : layerPositions2;
 994
 995       makeMissingMesh( linSide1, nodes, lp, myHelper );
 996       linSide1 = StdMeshers_FaceSide::New( nodes );
 997       makeMissingMesh( linSide2, nodes, lp, myHelper );
 998       linSide2 = StdMeshers_FaceSide::New( nodes );
 999     }
1000
1001     const UVPtStructVec& nodes2 = linSide2->GetUVPtStruct();
1002     centerNode = nodes2.back().node;
1003     centerUV   = nodes2.back().UV();
1004   }
1005
1006   list< TopoDS_Edge >::iterator ee = emptyEdges.begin();
1007   for ( ; ee != emptyEdges.end(); ++ee )
1008     TEdgeMarker::markEdge( *ee, aMesh.GetSubMesh( F ));
1009
1010   circSide->GetUVPtStruct(); // let sides take into account just computed nodes
1011   linSide1->GetUVPtStruct();
1012   linSide2->GetUVPtStruct();
1013
1014   FaceQuadStruct::Ptr quad( new FaceQuadStruct );
1015   quad->side.resize( 4 );
1016   quad->side[0] = circSide;
1017   quad->side[1] = linSide1;
1018   quad->side[2] = StdMeshers_FaceSide::New( circSide.get(), centerNode, &centerUV );
1019   quad->side[3] = linSide2;
1020   quad->face    = F;
1021
1022   myQuadList.push_back( quad );
1023
1024   // create quadrangles
1025   bool ok;
1026   if ( linSide1->NbPoints() == linSide2->NbPoints() )
1027     ok = StdMeshers_Quadrangle_2D::computeQuadDominant( aMesh, F, quad );
1028   else
1029     ok = StdMeshers_Quadrangle_2D::computeTriangles( aMesh, F, quad );
1030
1031   if ( helper.HasDegeneratedEdges() )
1032   {
1033     StdMeshers_Quadrangle_2D::myNeedSmooth = true;
1034     StdMeshers_Quadrangle_2D::smooth( quad );
1035   }
1036
1037   StdMeshers_Quadrangle_2D::myHelper = 0;
1038
1039   return ok;
1040 }
1041
1042 //================================================================================
1043 /*!
1044  * \brief Compute nodes on the radial edge
1045  * \retval int - nb of segments on the linSide
1046  */
1047 //================================================================================
1048
1049 int StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(StdMeshers_FaceSidePtr linSide,
1050                                                             vector< double > &     positions,
1051                                                             int*                   nbMeshedEdges,
1052                                                             bool                   useHalf)
1053 {
1054   // First, try to compute positions of layers
1055
1056   positions.clear();
1057
1058   SMESH_Mesh * mesh = myHelper->GetMesh();
1059
1060   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
1061   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
1062
1063   if ( !hyp1D && !nbLayers )
1064   {
1065     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
1066     // find a hyp usable by TNodeDistributor
1067     TopoDS_Shape edge = linSide->Edge(0);
1068     const SMESH_HypoFilter* hypKind =
1069       TNodeDistributor::GetDistributor(*mesh)->GetCompatibleHypoFilter(/*ignoreAux=*/true);
1070     hyp1D = mesh->GetHypothesis( edge, *hypKind, /*fromAncestors=*/true);
1071   }
1072   if ( hyp1D ) // try to compute with hyp1D
1073   {
1074     BRepAdaptor_CompCurve* curve = linSide->GetCurve3d();
1075     SMESHUtils::Deleter< BRepAdaptor_CompCurve > delCurve( curve );
1076     double f = curve->FirstParameter();
1077     double l = curve->LastParameter();
1078
1079     if ( useHalf )
1080       l = 0.5 * ( f + l ); // first part of linSide is used
1081
1082     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( positions, linSide->Edge(0),
1083                                                             *curve, f, l, *mesh, hyp1D ))
1084     {
1085       if ( myDistributionHypo ) { // bad hyp assigned
1086         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1087       }
1088       else {
1089         // bad hyp found, its Ok, lets try with default nb of segments
1090       }
1091     }
1092   }
1093
1094   if ( positions.empty() ) // try to use nb of layers
1095   {
1096     if ( !nbLayers )
1097       nbLayers = _gen->GetDefaultNbSegments();
1098
1099     if ( nbLayers )
1100     {
1101       positions.resize( nbLayers + 1 );
1102       for ( int z = 0; z < nbLayers; ++z )
1103         positions[ z ] = double( z )/ double( nbLayers );
1104       positions.back() = 1;
1105     }
1106   }
1107
1108   // Second, check presence of a mesh built by other algo on linSide
1109
1110   int nbEdgesComputed = 0;
1111   for ( int i = 0; i < linSide->NbEdges(); ++i )
1112   {
1113     nbEdgesComputed += ( !mesh->GetSubMesh( linSide->Edge(i))->IsEmpty() );
1114   }
1115
1116   if ( nbEdgesComputed == linSide->NbEdges() )
1117   {
1118     const UVPtStructVec& points = linSide->GetUVPtStruct();
1119     if ( points.size() >= 2 )
1120     {
1121       positions.resize( points.size() );
1122       for ( size_t i = 0; i < points.size(); ++i )
1123         positions[ i ] = points[i].normParam;
1124     }
1125   }
1126
1127   if ( nbMeshedEdges ) *nbMeshedEdges = nbEdgesComputed;
1128
1129   return positions.size();
1130 }
1131
1132
1133 //=======================================================================
1134 //function : Evaluate
1135 //purpose  :
1136 //=======================================================================
1137
1138 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh&         aMesh,
1139                                                 const TopoDS_Shape& aShape,
1140                                                 MapShapeNbElems&    aResMap)
1141 {
1142   if( aShape.ShapeType() != TopAbs_FACE ) {
1143     return false;
1144   }
1145   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1146   if( aResMap.count(sm) )
1147     return false;
1148
1149   vector<int>& aResVec =
1150     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1151
1152   myHelper = new SMESH_MesherHelper( aMesh );
1153   myHelper->SetSubShape( aShape );
1154   SMESHUtils::Deleter<SMESH_MesherHelper> helperDeleter( myHelper );
1155
1156   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1157
1158   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1159   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
1160   if( nbSides > 3 || nbSides < 1 )
1161     return false;
1162
1163   if ( algo1d->EvaluateCircularEdge( aMesh, circSide, aResMap ))
1164     return false;
1165
1166   vector< double > layerPositions; // [0,1]
1167
1168   // ------------------------------------------------------------------------------------------
1169   if ( nbSides == 1 )
1170   {
1171     const TopoDS_Face& F = TopoDS::Face( aShape );
1172
1173     const SMDS_MeshNode* circNode;
1174     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
1175
1176     StdMeshers_FaceSidePtr tmpSide =
1177       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
1178
1179     if ( !computeLayerPositions( tmpSide, layerPositions ))
1180       return false;
1181   }
1182   // ------------------------------------------------------------------------------------------
1183   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
1184   {
1185     if ( !computeLayerPositions( linSide1, layerPositions ))
1186       return false;
1187   }
1188   // ------------------------------------------------------------------------------------------
1189   else if ( nbSides == 2 )
1190   {
1191     // find positions of layers for the first half of linSide1
1192     if ( !computeLayerPositions( linSide1, layerPositions, 0, /*useHalf=*/true ))
1193       return false;
1194   }
1195   // ------------------------------------------------------------------------------------------
1196   else // nbSides == 3
1197   {
1198     if ( !computeLayerPositions(( linSide1->Length() > linSide2->Length() ) ? linSide1 : linSide2,
1199                                  layerPositions ))
1200       return false;
1201   }
1202
1203   bool isQuadratic = false;
1204   for ( TopExp_Explorer edge( aShape, TopAbs_EDGE ); edge.More() &&  !isQuadratic ; edge.Next() )
1205   {
1206     sm = aMesh.GetSubMesh( edge.Current() );
1207     vector<int>& nbElems = aResMap[ sm ];
1208     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1209       isQuadratic = nbElems[ SMDSEntity_Quad_Edge ];
1210   }
1211
1212   int nbCircSegments = 0;
1213   for ( int iE = 0; iE < circSide->NbEdges(); ++iE )
1214   {
1215     sm = aMesh.GetSubMesh( circSide->Edge( iE ));
1216     vector<int>& nbElems = aResMap[ sm ];
1217     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1218       nbCircSegments += ( nbElems[ SMDSEntity_Edge ] + nbElems[ SMDSEntity_Quad_Edge ]);
1219   }
1220
1221   int nbQuads = nbCircSegments * ( layerPositions.size() - 1 );
1222   int nbTria  = nbCircSegments;
1223   int nbNodes = ( nbCircSegments - 1 ) * ( layerPositions.size() - 2 );
1224   if ( isQuadratic )
1225   {
1226     nbNodes += (( nbCircSegments - 1 ) * ( layerPositions.size() - 1 ) + // radial
1227                 ( nbCircSegments     ) * ( layerPositions.size() - 2 )); // circular
1228     aResVec[SMDSEntity_Quad_Triangle  ] = nbTria;
1229     aResVec[SMDSEntity_Quad_Quadrangle] = nbQuads;
1230   }
1231   else
1232   {
1233     aResVec[SMDSEntity_Triangle  ] = nbTria;
1234     aResVec[SMDSEntity_Quadrangle] = nbQuads;
1235   }
1236   aResVec[SMDSEntity_Node] = nbNodes;
1237
1238   if ( linSide1 )
1239   {
1240     // evaluation for linSides
1241     vector<int> aResVec(SMDSEntity_Last, 0);
1242     if ( isQuadratic ) {
1243       aResVec[SMDSEntity_Node     ] = 2 * ( layerPositions.size() - 1 ) + 1;
1244       aResVec[SMDSEntity_Quad_Edge] = layerPositions.size() - 1;
1245     }
1246     else {
1247       aResVec[SMDSEntity_Node] = layerPositions.size() - 2;
1248       aResVec[SMDSEntity_Edge] = layerPositions.size() - 1;
1249     }
1250     sm = aMesh.GetSubMesh( linSide1->Edge(0) );
1251     aResMap[sm] = aResVec;
1252     if ( linSide2 )
1253     {
1254       sm = aMesh.GetSubMesh( linSide2->Edge(0) );
1255       aResMap[sm] = aResVec;
1256     }
1257   }
1258
1259   return true;
1260 }
1261
1262 //================================================================================
1263 /*!
1264  * \brief Return true if the algorithm can compute mesh on this shape
1265  */
1266 //================================================================================
1267
1268 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1269 {
1270   int nbFoundFaces = 0;
1271   for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
1272   {
1273     StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1274     int nbSides = analyseFace( exp.Current(), NULL, circSide, linSide1, linSide2, NULL );
1275     bool ok = ( 0 < nbSides && nbSides <= 3 &&
1276                 isCornerInsideCircle( circSide, linSide1, linSide2 ));
1277     if( toCheckAll  && !ok ) return false;
1278     if( !toCheckAll && ok  ) return true;
1279   }
1280   if( toCheckAll && nbFoundFaces != 0 ) return true;
1281   return false;
1282 }