src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

   1 // Copyright (C) 2007-2019  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     list<TopoDS_Edge>::iterator edge = edges.begin();
 230     while ( edge != edges.end() )
 231       if ( SMESH_Algo::isDegenerated( *edge ))
 232         edge = edges.erase( edge );
 233       else
 234         ++edge;
 235     int nbEdges = edges.size();
 236
 237     // find VERTEXes between continues EDGEs
 238     TopTools_MapOfShape contVV;
 239     if ( nbEdges > 1 )
 240     {
 241       TopoDS_Edge ePrev = edges.back();
 242       for ( edge = edges.begin(); edge != edges.end(); ++edge )
 243       {
 244         if ( SMESH_Algo::IsContinuous( ePrev, *edge ))
 245           contVV.Add( SMESH_MesherHelper::IthVertex( 0, *edge ));
 246         ePrev = *edge;
 247       }
 248     }
 249     // make edges start from a non-continues VERTEX
 250     if ( 1 < contVV.Extent() && contVV.Extent() < nbEdges )
 251     {
 252       while ( contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
 253         edges.splice( edges.end(), edges, edges.begin() );
 254     }
 255
 256     // make face sides
 257     TSideVector sides;
 258     while ( !edges.empty() )
 259     {
 260       list< TopoDS_Edge > sideEdges;
 261       sideEdges.splice( sideEdges.end(), edges, edges.begin() );
 262       while ( !edges.empty() &&
 263               contVV.Contains( SMESH_MesherHelper::IthVertex( 0, edges.front() )))
 264         sideEdges.splice( sideEdges.end(), edges, edges.begin() );
 265
 266       StdMeshers_FaceSidePtr side;
 267       if ( aMesh )
 268         side = StdMeshers_FaceSide::New( face, sideEdges, aMesh,
 269                                          /*isFwd=*/true, /*skipMedium=*/ true, helper );
 270       sides.push_back( side );
 271     }
 272
 273     if ( !aMesh ) // call from IsApplicable()
 274       return sides.size();
 275
 276     if ( sides.size() > 3 )
 277       return sides.size();
 278
 279     if ( nbWire == 2 && (( sides.size() != 2 ) ||
 280                          ( sides[0]->IsClosed() && sides[1]->IsClosed() ) ||
 281                          ( !sides[0]->IsClosed() && !sides[1]->IsClosed() )))
 282       return -1;
 283
 284     // detect an elliptic side
 285
 286     if ( sides.size() == 1 )
 287     {
 288       aCircSide = sides[0];
 289       return sides.size();
 290     }
 291
 292     // sort sides by deviation from a straight line
 293     multimap< double, int > deviation2sideInd;
 294     const double nbSamples = 7;
 295     for ( size_t iS = 0; iS < sides.size(); ++iS )
 296     {
 297       gp_Pnt pf = BRep_Tool::Pnt( sides[iS]->FirstVertex() );
 298       gp_Pnt pl = BRep_Tool::Pnt( sides[iS]->LastVertex() );
 299       gp_Vec v1( pf, pl );
 300       double v1Len = v1.Magnitude();
 301       if ( v1Len < std::numeric_limits< double >::min() )
 302       {
 303         deviation2sideInd.insert( make_pair( sides[iS]->Length(), iS )); // the side seems closed
 304         continue;
 305       }
 306       double devia = 0;
 307       for ( int i = 0; i < nbSamples; ++i )
 308       {
 309         gp_Pnt pi( sides[iS]->Value3d(( i + 1 ) / nbSamples ));
 310         gp_Vec vi( pf, pi );
 311         double h = 0.5 * v1.Crossed( vi ).Magnitude() / v1Len;
 312         devia = Max( devia, h );
 313       }
 314       deviation2sideInd.insert( make_pair( devia, iS ));
 315     }
 316     double maxDevi = deviation2sideInd.rbegin()->first;
 317     if ( maxDevi < 1e-7 && sides.size() == 3 )
 318     {
 319       // a triangle FACE; use a side with the most outstanding length as an elliptic one
 320       deviation2sideInd.clear();
 321       multimap< double, int > len2sideInd;
 322       for ( size_t iS = 0; iS < sides.size(); ++iS )
 323         len2sideInd.insert( make_pair( sides[iS]->Length(), iS ));
 324
 325       multimap< double, int >::iterator l2i = len2sideInd.begin();
 326       double len0 = l2i->first;
 327       double len1 = (++l2i)->first;
 328       double len2 = (++l2i)->first;
 329       if ( len1 - len0 > len2 - len1 )
 330         deviation2sideInd.insert( make_pair( 0., len2sideInd.begin()->second ));
 331       else
 332         deviation2sideInd.insert( make_pair( 0., len2sideInd.rbegin()->second ));
 333     }
 334
 335     int iCirc = deviation2sideInd.rbegin()->second;
 336     aCircSide = sides[ iCirc ];
 337     aLinSide1 = sides[( iCirc + 1 ) % sides.size() ];
 338     if ( sides.size() > 2 )
 339     {
 340       aLinSide2 = sides[( iCirc + 2 ) % sides.size() ];
 341       aLinSide2->Reverse(); // to be "parallel" to aLinSide1
 342     }
 343
 344     if (( nbWire == 2 && aLinSide1 ) &&
 345         ( aLinSide1->Edge(0).Orientation() == TopAbs_INTERNAL ) &&
 346         ( aCircSide->IsClosed() ))
 347     {
 348       // assure that aCircSide starts at aLinSide1
 349       TopoDS_Vertex v0 = aLinSide1->FirstVertex();
 350       TopoDS_Vertex v1 = aLinSide1->LastVertex();
 351       if ( ! aCircSide->FirstVertex().IsSame( v0 ) &&
 352            ! aCircSide->FirstVertex().IsSame( v1 ))
 353       {
 354         int iE = 0;
 355         for ( ; iE < aCircSide->NbEdges(); ++iE )
 356           if ( aCircSide->FirstVertex(iE).IsSame( v0 ) ||
 357                aCircSide->FirstVertex(iE).IsSame( v1 ))
 358             break;
 359         if ( iE == aCircSide->NbEdges() )
 360           return -2;
 361
 362         edges.clear();
 363         for ( int i = 0; i < aCircSide->NbEdges(); ++i, ++iE )
 364           edges.push_back( aCircSide->Edge( iE % aCircSide->NbEdges() ));
 365
 366         aCircSide = StdMeshers_FaceSide::New( face, edges, aMesh,
 367                                               /*isFwd=*/true, /*skipMedium=*/ true, helper );
 368       }
 369     }
 370
 371     return sides.size();
 372   }
 373
 374   //================================================================================
 375   /*!
 376    * \brief Checks if the common vertex between LinSide's lies inside the circle
 377    *  and not outside
 378    *  \return bool - false if there are 3 EDGEs and the corner is outside
 379    */
 380   //================================================================================
 381
 382   bool isCornerInsideCircle(const StdMeshers_FaceSidePtr& CircSide,
 383                             const StdMeshers_FaceSidePtr& LinSide1,
 384                             const StdMeshers_FaceSidePtr& LinSide2)
 385   {
 386     // if ( CircSide && LinSide1 && LinSide2 )
 387     // {
 388     //   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircSide ));
 389     //   TopoDS_Vertex aCommonV;
 390     //   if ( !aCirc.IsNull() &&
 391     //        TopExp::CommonVertex( LinSide1, LinSide2, aCommonV ))
 392     //   {
 393     //     gp_Pnt aCommonP = BRep_Tool::Pnt( aCommonV );
 394     //     gp_Pnt  aCenter = aCirc->Location();
 395     //     double     dist = aCenter.Distance( aCommonP );
 396     //     return dist < 0.1 * aCirc->Radius();
 397     //   }
 398     // }
 399     return true;
 400   }
 401
 402   //================================================================================
 403   /*!
 404    * \brief Create an EDGE connecting the ellipse center with the most distant point
 405    *        of the ellipse.
 406    *  \param [in] circSide - the elliptic side
 407   *  \param [in] face - the FACE
 408   *  \param [out] circNode - a node on circSide most distant from the center
 409   *  \return TopoDS_Edge - the create EDGE
 410   */
 411   //================================================================================
 412
 413   TopoDS_Edge makeEdgeToCenter( StdMeshers_FaceSidePtr& circSide,
 414                                 const TopoDS_Face&      face,
 415                                 const SMDS_MeshNode*&   circNode)
 416   {
 417     // find the center and a point most distant from it
 418
 419     double maxDist = 0, normPar;
 420     gp_XY uv1, uv2;
 421     for ( int i = 0; i < 32; ++i )
 422     {
 423       double    u = 0.5 * i / 32.;
 424       gp_Pnt2d p1 = circSide->Value2d( u );
 425       gp_Pnt2d p2 = circSide->Value2d( u + 0.5 );
 426       double dist = p1.SquareDistance( p2 );
 427       if ( dist > maxDist )
 428       {
 429         maxDist = dist;
 430         uv1 = p1.XY();
 431         uv2 = p2.XY();
 432         normPar = u;
 433       }
 434     }
 435     gp_XY center = 0.5 * ( uv1 + uv2 );
 436     double   len = 0.5 * Sqrt( maxDist );
 437     bool  isCirc = ( Abs( len - circSide->Value2d( 0 ).Distance( center )) < 1e-3 * len );
 438
 439     // find a node closest to the most distant point
 440
 441     size_t iDist = 0;
 442     const UVPtStructVec& circNodes = circSide->GetUVPtStruct();
 443     if ( !isCirc )
 444     {
 445       double minDist = 1e100;
 446       for ( size_t i = 0; i <= circNodes.size(); ++i )
 447       {
 448         double dist = Abs( circNodes[i].normParam - normPar );
 449         if ( dist < minDist )
 450         {
 451           iDist = i;
 452           minDist = dist;
 453         }
 454       }
 455     }
 456     circNode = circNodes[iDist].node;
 457     uv1      = circNodes[iDist].UV();
 458     len      = ( uv1 - center ).Modulus();
 459
 460     // make the EDGE between the most distant point and the center
 461
 462     Handle(Geom2d_Line) line = new Geom2d_Line( uv1, gp_Dir2d( center - uv1 ) );
 463     Handle(Geom2d_Curve) pcu = new Geom2d_TrimmedCurve( line, 0, len );
 464
 465     Handle(Geom_Surface) surface = BRep_Tool::Surface( face );
 466     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pcu, surface, 0, len );
 467
 468     BRep_Builder().UpdateEdge( edge, pcu, face, 1e-7 );
 469     ShapeFix_Edge().FixAddCurve3d( edge );
 470
 471
 472     // assure that circSide starts at circNode
 473     if ( iDist != 0  &&  iDist != circNodes.size()-1 )
 474     {
 475       // create a new circSide
 476       UVPtStructVec nodesNew;
 477       nodesNew.reserve( circNodes.size() );
 478       nodesNew.insert( nodesNew.end(), circNodes.begin() + iDist, circNodes.end() );
 479       nodesNew.insert( nodesNew.end(), circNodes.begin() + 1, circNodes.begin() + iDist + 1 );
 480       circSide = StdMeshers_FaceSide::New( nodesNew );
 481     }
 482
 483     return edge;
 484   }
 485
 486   //================================================================================
 487   /*!
 488    * \brief Set nodes existing on a linSide to UVPtStructVec and create missing nodes
 489    *        corresponding to layerPositions
 490    */
 491   //================================================================================
 492
 493   void makeMissingMesh( StdMeshers_FaceSidePtr& linSide,
 494                         UVPtStructVec&          nodes,
 495                         const vector< double >& layerPositions,
 496                         SMESH_MesherHelper*     helper )
 497   {
 498     // tolerance to compare normParam
 499     double tol = 1e100;
 500     for ( size_t i = 1; i < layerPositions.size(); ++i )
 501       tol = Min( tol, layerPositions[i] - layerPositions[i-1] );
 502     tol *= 0.05;
 503
 504     // merge existing nodes with layerPositions into UVPtStructVec
 505     // ------------------------------------------------------------
 506
 507     const UVPtStructVec& exiNodes = linSide->GetUVPtStruct();
 508     nodes.clear();
 509     nodes.reserve( layerPositions.size() + exiNodes.size() );
 510     vector< double >::const_iterator pos = layerPositions.begin(), posEnd = layerPositions.end();
 511     for ( size_t i = 0; i < exiNodes.size(); ++i )
 512     {
 513       switch ( exiNodes[i].node->GetPosition()->GetTypeOfPosition() )
 514       {
 515       case SMDS_TOP_VERTEX:
 516       {
 517         // allocate UVPtStruct's for non-existing nodes
 518         while ( pos != posEnd && *pos < exiNodes[i].normParam - tol )
 519         {
 520           UVPtStruct uvPS;
 521           uvPS.normParam = *pos++;
 522           nodes.push_back( uvPS );
 523         }
 524         // save existing node on a VERTEX
 525         nodes.push_back( exiNodes[i] );
 526         break;
 527       }
 528       case SMDS_TOP_EDGE:
 529       {
 530         // save existing nodes on an EDGE
 531         while ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 1 )
 532         {
 533           nodes.push_back( exiNodes[i++] );
 534         }
 535         // save existing node on a VERTEX
 536         if ( i < exiNodes.size() && exiNodes[i].node->GetPosition()->GetDim() == 0 )
 537         {
 538           nodes.push_back( exiNodes[i] );
 539         }
 540         break;
 541       }
 542       default:;
 543       }
 544
 545       // skip layer positions covered by saved nodes
 546       while ( pos != posEnd && *pos < nodes.back().normParam + tol )
 547       {
 548         ++pos;
 549       }
 550     }
 551     // allocate UVPtStruct's for the rest non-existing nodes
 552     while ( pos != posEnd )
 553     {
 554       UVPtStruct uvPS;
 555       uvPS.normParam = *pos++;
 556       nodes.push_back( uvPS );
 557     }
 558
 559     // create missing nodes
 560     // ---------------------
 561
 562     SMESHDS_Mesh * meshDS = helper->GetMeshDS();
 563     const TopoDS_Face& F  = TopoDS::Face( helper->GetSubShape() );
 564     Handle(ShapeAnalysis_Surface) surface = helper->GetSurface( F );
 565
 566     helper->SetElementsOnShape( false ); // we create nodes on EDGEs, not on the FACE
 567
 568     for ( size_t i = 0; i < nodes.size(); ++i )
 569     {
 570       if ( nodes[ i ].node ) continue;
 571
 572       gp_Pnt2d uv = linSide->Value2d( nodes[i].normParam );
 573       gp_Pnt  xyz = surface->Value( uv.X(), uv.Y() );
 574
 575       nodes[ i ].SetUV( uv.XY() );
 576       nodes[ i ].node = helper->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
 577     }
 578
 579     // set nodes on VERTEXes
 580     for ( int iE = 0; iE < linSide->NbEdges(); ++iE )
 581     {
 582       TopoDS_Vertex v = linSide->LastVertex( iE );
 583       if ( SMESH_Algo::VertexNode( v, meshDS ))
 584         continue;
 585
 586       double normPar = linSide->LastParameter( iE );
 587       size_t i = 0;
 588       while ( nodes[ i ].normParam < normPar )
 589         ++i;
 590       if (( nodes[ i ].normParam - normPar ) > ( normPar - nodes[ i-1 ].normParam ))
 591         --i;
 592       meshDS->SetNodeOnVertex( nodes[ i ].node, v );
 593     }
 594
 595     // set nodes on EDGEs
 596     int edgeID;
 597     for ( size_t i = 0; i < nodes.size(); ++i )
 598     {
 599       if ( nodes[ i ].node->getshapeId() > 0 ) continue;
 600
 601       double u = linSide->Parameter( nodes[i].normParam, edgeID );
 602       meshDS->SetNodeOnEdge( nodes[ i ].node, edgeID, u );
 603     }
 604
 605     // create segments
 606     for ( size_t i = 1; i < nodes.size(); ++i )
 607     {
 608       if ( meshDS->FindEdge( nodes[i].node, nodes[i-1].node )) continue;
 609
 610       const SMDS_MeshElement* seg = helper->AddEdge( nodes[i].node, nodes[i-1].node );
 611
 612       double normParam = 0.5 * ( nodes[i].normParam + nodes[i-1].normParam );
 613       edgeID = linSide->EdgeID( linSide->EdgeIndex( normParam ));
 614       meshDS->SetMeshElementOnShape( seg, edgeID );
 615     }
 616
 617     helper->SetElementsOnShape( true );
 618
 619   } // end makeMissingMesh()
 620
 621 //================================================================================
 622 //================================================================================
 623 /*!
 624  * \brief Class computing layers distribution using data of
 625  *        StdMeshers_LayerDistribution hypothesis
 626  */
 627 //================================================================================
 628 //================================================================================
 629
 630 class TNodeDistributor: public StdMeshers_Regular_1D
 631 {
 632   list <const SMESHDS_Hypothesis *> myUsedHyps;
 633 public:
 634   // -----------------------------------------------------------------------------
 635   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 636   {
 637     const int myID = -1001;
 638     TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
 639     if ( !myHyp )
 640       myHyp = new TNodeDistributor( myID, aMesh.GetGen() );
 641     return myHyp;
 642   }
 643   // -----------------------------------------------------------------------------
 644   //! Computes distribution of nodes on a straight line ending at pIn and pOut
 645   bool Compute( vector< double > &      positions,
 646                 const TopoDS_Edge&      edge,
 647                 Adaptor3d_Curve&        curve,
 648                 double                  f,
 649                 double                  l,
 650                 SMESH_Mesh&             mesh,
 651                 const SMESH_Hypothesis* hyp1d)
 652   {
 653     if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
 654
 655     myUsedHyps.clear();
 656     myUsedHyps.push_back( hyp1d );
 657
 658     SMESH_Hypothesis::Hypothesis_Status aStatus;
 659     if ( !StdMeshers_Regular_1D::CheckHypothesis( mesh, edge, aStatus ))
 660       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 661                     "with LayerDistribution hypothesis");
 662
 663     double len = GCPnts_AbscissaPoint::Length( curve, f, l );
 664
 665     list< double > params;
 666     if ( !StdMeshers_Regular_1D::computeInternalParameters( mesh, curve, len, f, l, params, false ))
 667       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 668
 669     params.push_front( f );
 670     params.push_back ( l );
 671     positions.clear();
 672     positions.reserve( params.size() );
 673     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 674       positions.push_back(( *itU - f ) / ( l - f ));
 675     return true;
 676   }
 677   // -----------------------------------------------------------------------------
 678   //! Make mesh on an edge using assigned 1d hyp or default nb of segments
 679   bool ComputeCircularEdge( SMESH_Mesh&                   aMesh,
 680                             const StdMeshers_FaceSidePtr& aSide )
 681   {
 682     bool ok = true;
 683     for ( int i = 0; i < aSide->NbEdges(); ++i )
 684     {
 685       const TopoDS_Edge& edge = aSide->Edge( i );
 686       _gen->Compute( aMesh, edge );
 687       SMESH_subMesh *sm = aMesh.GetSubMesh( edge );
 688       if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
 689       {
 690         // find any 1d hyp assigned (there can be a hyp w/o algo)
 691         myUsedHyps = SMESH_Algo::GetUsedHypothesis( aMesh, edge, /*ignoreAux=*/true );
 692         Hypothesis_Status aStatus;
 693         if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 694         {
 695           // no valid 1d hyp assigned, use default nb of segments
 696           _hypType                    = NB_SEGMENTS;
 697           _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 698           _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 699         }
 700         ok &= StdMeshers_Regular_1D::Compute( aMesh, edge );
 701       }
 702     }
 703     return ok;
 704   }
 705   // -----------------------------------------------------------------------------
 706   //! Make mesh on an edge using assigned 1d hyp or default nb of segments
 707   bool EvaluateCircularEdge(SMESH_Mesh&                  aMesh,
 708                             const StdMeshers_FaceSidePtr aSide,
 709                             MapShapeNbElems&             aResMap)
 710   {
 711     bool ok = true;
 712     for ( int i = 0; i < aSide->NbEdges(); ++i )
 713     {
 714       const TopoDS_Edge& anEdge = aSide->Edge( i );
 715       _gen->Evaluate( aMesh, anEdge, aResMap );
 716       if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
 717         continue;
 718
 719       // find any 1d hyp assigned
 720       myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
 721       Hypothesis_Status aStatus;
 722       if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
 723       {
 724         // no valid 1d hyp assigned, use default nb of segments
 725         _hypType                    = NB_SEGMENTS;
 726         _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 727         _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 728       }
 729       ok &= StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
 730     }
 731     return ok;
 732   }
 733 protected:
 734   // -----------------------------------------------------------------------------
 735   TNodeDistributor( int hypId, SMESH_Gen* gen)
 736     : StdMeshers_Regular_1D( hypId, gen)
 737   {
 738   }
 739   // -----------------------------------------------------------------------------
 740   virtual const list <const SMESHDS_Hypothesis *> &
 741     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 742   {
 743     return myUsedHyps;
 744   }
 745   // -----------------------------------------------------------------------------
 746 };
 747 }
 748
 749 //=======================================================================
 750 /*!
 751  * \brief Allow algo to do something after persistent restoration
 752  * \param subMesh - restored submesh
 753  *
 754  * call TEdgeMarker::markEdge()
 755  */
 756 //=======================================================================
 757
 758 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
 759 {
 760   if ( !faceSubMesh->IsEmpty() )
 761     SetEventListener( faceSubMesh );
 762 }
 763
 764 //=======================================================================
 765 /*!
 766  * \brief Sets event listener to a submesh
 767  * \param subMesh - submesh where algo is set
 768  *
 769  * This method is called when a submesh gets HYP_OK algo_state.
 770  */
 771 //=======================================================================
 772
 773 void StdMeshers_RadialQuadrangle_1D2D::SetEventListener(SMESH_subMesh* faceSubMesh)
 774 {
 775   for ( TopExp_Explorer e( faceSubMesh->GetSubShape(), TopAbs_EDGE ); e.More(); e.Next() )
 776   {
 777     TEdgeMarker::markEdge( TopoDS::Edge( e.Current() ), faceSubMesh );
 778   }
 779 }
 780
 781 //=======================================================================
 782 //function : Compute
 783 //purpose  :
 784 //=======================================================================
 785
 786 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 787                                                const TopoDS_Shape& aShape)
 788 {
 789   SMESH_MesherHelper helper( aMesh );
 790   StdMeshers_Quadrangle_2D::myHelper     = & helper;
 791   StdMeshers_Quadrangle_2D::myNeedSmooth = false;
 792   StdMeshers_Quadrangle_2D::myCheckOri   = false;
 793   StdMeshers_Quadrangle_2D::myQuadList.clear();
 794
 795   myHelper->SetSubShape( aShape );
 796   myHelper->SetElementsOnShape( true );
 797
 798   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
 799   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
 800   if( nbSides > 3 || nbSides < 1 )
 801     return error("The face must be a full ellipse or a part of ellipse (i.e. the number "
 802                  "of edges is less or equal to 3 and one of them is an ellipse curve)");
 803
 804   // get not yet computed EDGEs
 805   list< TopoDS_Edge > emptyEdges;
 806   for ( TopExp_Explorer e( aShape, TopAbs_EDGE ); e.More(); e.Next() )
 807   {
 808     if ( aMesh.GetSubMesh( e.Current() )->IsEmpty() )
 809       emptyEdges.push_back( TopoDS::Edge( e.Current() ));
 810   }
 811
 812   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
 813
 814   if ( !algo1d->ComputeCircularEdge( aMesh, circSide ))
 815     return error( algo1d->GetComputeError() );
 816
 817
 818   TopoDS_Face F = TopoDS::Face(aShape);
 819   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 820
 821   myHelper->IsQuadraticSubMesh( aShape );
 822
 823   vector< double > layerPositions; // [0,1]
 824
 825   const SMDS_MeshNode* centerNode = 0;
 826   gp_Pnt2d             centerUV(0,0);
 827
 828   // ------------------------------------------------------------------------------------------
 829   if ( nbSides == 1 ) // full ellipse
 830   {
 831     const SMDS_MeshNode* circNode;
 832     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
 833
 834     StdMeshers_FaceSidePtr tmpSide =
 835       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
 836
 837     if ( !computeLayerPositions( tmpSide, layerPositions ))
 838       return false;
 839
 840     UVPtStructVec nodes( layerPositions.size() );
 841     nodes[0].node = circNode;
 842     for ( size_t i = 0; i < layerPositions.size(); ++i )
 843     {
 844       gp_Pnt2d uv = tmpSide->Value2d( layerPositions[i] );
 845       gp_Pnt  xyz = S->Value( uv.X(), uv.Y() );
 846
 847       nodes[ i ].SetUV( uv.XY() );
 848       nodes[ i ].normParam = layerPositions[i];
 849       if ( i )
 850         nodes[ i ].node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, uv.X(), uv.Y() );
 851     }
 852
 853     linSide1 = StdMeshers_FaceSide::New( nodes );
 854     linSide2 = StdMeshers_FaceSide::New( nodes );
 855
 856     centerNode = nodes.back().node;
 857     centerUV   = nodes.back().UV();
 858   }
 859   // ------------------------------------------------------------------------------------------
 860   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
 861   {
 862     // full ellipse with an internal radial side
 863
 864     // eliminate INTERNAL orientation
 865     list< TopoDS_Edge > edges;
 866     for ( int iE = 0; iE < linSide1->NbEdges(); ++iE )
 867     {
 868       edges.push_back( linSide1->Edge( iE ));
 869       edges.back().Orientation( TopAbs_FORWARD );
 870     }
 871
 872     // orient the internal side
 873     bool isVIn0Shared = false;
 874     TopoDS_Vertex vIn0 = myHelper->IthVertex( 0, edges.front() );
 875     for ( int iE = 0; iE < circSide->NbEdges() && !isVIn0Shared; ++iE )
 876       isVIn0Shared = vIn0.IsSame( circSide->FirstVertex( iE ));
 877
 878     linSide1 = StdMeshers_FaceSide::New( F, edges, &aMesh,
 879                                          /*isFrw=*/isVIn0Shared, /*skipMedium=*/true, myHelper );
 880
 881     int nbMeshedEdges;
 882     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges ))
 883       return false;
 884
 885     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 886     UVPtStructVec nodes;
 887     if ( nbMeshedEdges != linSide1->NbEdges() )
 888       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 889     else
 890       nodes = linSide1->GetUVPtStruct();
 891
 892     linSide1 = StdMeshers_FaceSide::New( nodes );
 893     linSide2 = StdMeshers_FaceSide::New( nodes );
 894
 895     centerNode = nodes.back().node;
 896     centerUV   = nodes.back().UV();
 897   }
 898   // ------------------------------------------------------------------------------------------
 899   else if ( nbSides == 2 )
 900   {
 901     // find positions of layers for the first half of linSide1
 902     int nbMeshedEdges;
 903     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges, /*useHalf=*/true ))
 904       return false;
 905
 906     // make positions for the whole linSide1
 907     for ( size_t i = 0; i < layerPositions.size(); ++i )
 908     {
 909       layerPositions[i] *= 0.5;
 910     }
 911     layerPositions.reserve( layerPositions.size() * 2 );
 912     for ( int nb = layerPositions.size()-1; nb > 0; --nb )
 913       layerPositions.push_back( layerPositions.back() + layerPositions[nb] - layerPositions[nb-1] );
 914
 915     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 916     UVPtStructVec nodes;
 917     if ( nbMeshedEdges != linSide1->NbEdges() )
 918       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 919     else
 920       nodes = linSide1->GetUVPtStruct();
 921
 922     // find a central node
 923     size_t i = 0;
 924     while ( nodes[ i ].normParam < 0.5 ) ++i;
 925     if (( nodes[ i ].normParam - 0.5 ) > ( 0.5 - nodes[ i-1 ].normParam )) --i;
 926
 927     // distribute nodes between two linear sides
 928     UVPtStructVec nodes2( nodes.rbegin(), nodes.rbegin() + nodes.size() - i );
 929     nodes.resize( i + 1 );
 930
 931     linSide1 = StdMeshers_FaceSide::New( nodes );
 932     linSide2 = StdMeshers_FaceSide::New( nodes2 );
 933
 934     centerNode = nodes.back().node;
 935     centerUV   = nodes.back().UV();
 936   }
 937   // ------------------------------------------------------------------------------------------
 938   else // nbSides == 3
 939   {
 940     // one curve must be a part of ellipse and 2 other curves must be segments of line
 941
 942     int nbMeshedEdges1, nbMeshedEdges2;
 943     vector< double > layerPositions2;
 944     bool ok1 = computeLayerPositions( linSide1, layerPositions,  &nbMeshedEdges1 );
 945     bool ok2 = computeLayerPositions( linSide2, layerPositions2, &nbMeshedEdges2 );
 946     if ( !ok1 && !ok2 )
 947       return false;
 948
 949     bool linSide1Computed = ( nbMeshedEdges1 == linSide1->NbEdges() );
 950     bool linSide2Computed = ( nbMeshedEdges2 == linSide2->NbEdges() );
 951
 952     UVPtStructVec nodes;
 953
 954     if ( linSide1Computed && !linSide2Computed )
 955     {
 956       // use layer positions of linSide1 to mesh linSide2
 957       makeMissingMesh( linSide2, nodes, layerPositions, myHelper );
 958       linSide2 = StdMeshers_FaceSide::New( nodes );
 959     }
 960     else if ( linSide2Computed && !linSide1Computed )
 961     {
 962       // use layer positions of linSide2 to mesh linSide1
 963       makeMissingMesh( linSide1, nodes, layerPositions2, myHelper );
 964       linSide1 = StdMeshers_FaceSide::New( nodes );
 965     }
 966     else if ( !linSide2Computed && !linSide1Computed )
 967     {
 968       // use layer positions of a longer side to mesh the shorter side
 969       vector< double >& lp =
 970         ( linSide1->Length() > linSide2->Length() ) ? layerPositions : layerPositions2;
 971
 972       makeMissingMesh( linSide1, nodes, lp, myHelper );
 973       linSide1 = StdMeshers_FaceSide::New( nodes );
 974       makeMissingMesh( linSide2, nodes, lp, myHelper );
 975       linSide2 = StdMeshers_FaceSide::New( nodes );
 976     }
 977
 978     const UVPtStructVec& nodes2 = linSide2->GetUVPtStruct();
 979     centerNode = nodes2.back().node;
 980     centerUV   = nodes2.back().UV();
 981   }
 982
 983   list< TopoDS_Edge >::iterator ee = emptyEdges.begin();
 984   for ( ; ee != emptyEdges.end(); ++ee )
 985     TEdgeMarker::markEdge( *ee, aMesh.GetSubMesh( F ));
 986
 987   circSide->GetUVPtStruct(); // let sides take into account just computed nodes
 988   linSide1->GetUVPtStruct();
 989   linSide2->GetUVPtStruct();
 990
 991   FaceQuadStruct::Ptr quad( new FaceQuadStruct );
 992   quad->side.resize( 4 );
 993   quad->side[0] = circSide;
 994   quad->side[1] = linSide1;
 995   quad->side[2] = StdMeshers_FaceSide::New( circSide.get(), centerNode, &centerUV );
 996   quad->side[3] = linSide2;
 997
 998   myQuadList.push_back( quad );
 999
1000   // create quadrangles
1001   bool ok;
1002   if ( linSide1->NbPoints() == linSide2->NbPoints() )
1003     ok = StdMeshers_Quadrangle_2D::computeQuadDominant( aMesh, F, quad );
1004   else
1005     ok = StdMeshers_Quadrangle_2D::computeTriangles( aMesh, F, quad );
1006
1007   StdMeshers_Quadrangle_2D::myHelper = 0;
1008
1009   return ok;
1010 }
1011
1012 //================================================================================
1013 /*!
1014  * \brief Compute nodes on the radial edge
1015  * \retval int - nb of segments on the linSide
1016  */
1017 //================================================================================
1018
1019 int StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(StdMeshers_FaceSidePtr linSide,
1020                                                             vector< double > &     positions,
1021                                                             int*                   nbMeshedEdges,
1022                                                             bool                   useHalf)
1023 {
1024   // First, try to compute positions of layers
1025
1026   positions.clear();
1027
1028   SMESH_Mesh * mesh = myHelper->GetMesh();
1029
1030   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
1031   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
1032
1033   if ( !hyp1D && !nbLayers )
1034   {
1035     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
1036     // find a hyp usable by TNodeDistributor
1037     TopoDS_Shape edge = linSide->Edge(0);
1038     const SMESH_HypoFilter* hypKind =
1039       TNodeDistributor::GetDistributor(*mesh)->GetCompatibleHypoFilter(/*ignoreAux=*/true);
1040     hyp1D = mesh->GetHypothesis( edge, *hypKind, /*fromAncestors=*/true);
1041   }
1042   if ( hyp1D ) // try to compute with hyp1D
1043   {
1044     BRepAdaptor_CompCurve* curve = linSide->GetCurve3d();
1045     SMESHUtils::Deleter< BRepAdaptor_CompCurve > delCurve( curve );
1046     double f = curve->FirstParameter();
1047     double l = curve->LastParameter();
1048
1049     if ( useHalf )
1050       l = 0.5 * ( f + l ); // first part of linSide is used
1051
1052     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( positions, linSide->Edge(0),
1053                                                             *curve, f, l, *mesh, hyp1D ))
1054     {
1055       if ( myDistributionHypo ) { // bad hyp assigned
1056         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1057       }
1058       else {
1059         // bad hyp found, its Ok, lets try with default nb of segments
1060       }
1061     }
1062   }
1063
1064   if ( positions.empty() ) // try to use nb of layers
1065   {
1066     if ( !nbLayers )
1067       nbLayers = _gen->GetDefaultNbSegments();
1068
1069     if ( nbLayers )
1070     {
1071       positions.resize( nbLayers + 1 );
1072       for ( int z = 0; z < nbLayers; ++z )
1073         positions[ z ] = double( z )/ double( nbLayers );
1074       positions.back() = 1;
1075     }
1076   }
1077
1078   // Second, check presence of a mesh built by other algo on linSide
1079
1080   int nbEdgesComputed = 0;
1081   for ( int i = 0; i < linSide->NbEdges(); ++i )
1082   {
1083     nbEdgesComputed += ( !mesh->GetSubMesh( linSide->Edge(i))->IsEmpty() );
1084   }
1085
1086   if ( nbEdgesComputed == linSide->NbEdges() )
1087   {
1088     const UVPtStructVec& points = linSide->GetUVPtStruct();
1089     if ( points.size() >= 2 )
1090     {
1091       positions.resize( points.size() );
1092       for ( size_t i = 0; i < points.size(); ++i )
1093         positions[ i ] = points[i].normParam;
1094     }
1095   }
1096
1097   if ( nbMeshedEdges ) *nbMeshedEdges = nbEdgesComputed;
1098
1099   return positions.size();
1100 }
1101
1102
1103 //=======================================================================
1104 //function : Evaluate
1105 //purpose  :
1106 //=======================================================================
1107
1108 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh&         aMesh,
1109                                                 const TopoDS_Shape& aShape,
1110                                                 MapShapeNbElems&    aResMap)
1111 {
1112   if( aShape.ShapeType() != TopAbs_FACE ) {
1113     return false;
1114   }
1115   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1116   if( aResMap.count(sm) )
1117     return false;
1118
1119   vector<int>& aResVec =
1120     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1121
1122   myHelper = new SMESH_MesherHelper( aMesh );
1123   myHelper->SetSubShape( aShape );
1124   SMESHUtils::Deleter<SMESH_MesherHelper> helperDeleter( myHelper );
1125
1126   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1127
1128   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1129   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
1130   if( nbSides > 3 || nbSides < 1 )
1131     return false;
1132
1133   if ( algo1d->EvaluateCircularEdge( aMesh, circSide, aResMap ))
1134     return false;
1135
1136   vector< double > layerPositions; // [0,1]
1137
1138   // ------------------------------------------------------------------------------------------
1139   if ( nbSides == 1 )
1140   {
1141     const TopoDS_Face& F = TopoDS::Face( aShape );
1142
1143     const SMDS_MeshNode* circNode;
1144     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
1145
1146     StdMeshers_FaceSidePtr tmpSide =
1147       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
1148
1149     if ( !computeLayerPositions( tmpSide, layerPositions ))
1150       return false;
1151   }
1152   // ------------------------------------------------------------------------------------------
1153   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
1154   {
1155     if ( !computeLayerPositions( linSide1, layerPositions ))
1156       return false;
1157   }
1158   // ------------------------------------------------------------------------------------------
1159   else if ( nbSides == 2 )
1160   {
1161     // find positions of layers for the first half of linSide1
1162     if ( !computeLayerPositions( linSide1, layerPositions, 0, /*useHalf=*/true ))
1163       return false;
1164   }
1165   // ------------------------------------------------------------------------------------------
1166   else // nbSides == 3
1167   {
1168     if ( !computeLayerPositions(( linSide1->Length() > linSide2->Length() ) ? linSide1 : linSide2,
1169                                  layerPositions ))
1170       return false;
1171   }
1172
1173   bool isQuadratic = false;
1174   for ( TopExp_Explorer edge( aShape, TopAbs_EDGE ); edge.More() &&  !isQuadratic ; edge.Next() )
1175   {
1176     sm = aMesh.GetSubMesh( edge.Current() );
1177     vector<int>& nbElems = aResMap[ sm ];
1178     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1179       isQuadratic = nbElems[ SMDSEntity_Quad_Edge ];
1180   }
1181
1182   int nbCircSegments = 0;
1183   for ( int iE = 0; iE < circSide->NbEdges(); ++iE )
1184   {
1185     sm = aMesh.GetSubMesh( circSide->Edge( iE ));
1186     vector<int>& nbElems = aResMap[ sm ];
1187     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1188       nbCircSegments += ( nbElems[ SMDSEntity_Edge ] + nbElems[ SMDSEntity_Quad_Edge ]);
1189   }
1190
1191   int nbQuads = nbCircSegments * ( layerPositions.size() - 1 );
1192   int nbTria  = nbCircSegments;
1193   int nbNodes = ( nbCircSegments - 1 ) * ( layerPositions.size() - 2 );
1194   if ( isQuadratic )
1195   {
1196     nbNodes += (( nbCircSegments - 1 ) * ( layerPositions.size() - 1 ) + // radial
1197                 ( nbCircSegments     ) * ( layerPositions.size() - 2 )); // circular
1198     aResVec[SMDSEntity_Quad_Triangle  ] = nbTria;
1199     aResVec[SMDSEntity_Quad_Quadrangle] = nbQuads;
1200   }
1201   else
1202   {
1203     aResVec[SMDSEntity_Triangle  ] = nbTria;
1204     aResVec[SMDSEntity_Quadrangle] = nbQuads;
1205   }
1206   aResVec[SMDSEntity_Node] = nbNodes;
1207
1208   if ( linSide1 )
1209   {
1210     // evaluation for linSides
1211     vector<int> aResVec(SMDSEntity_Last, 0);
1212     if ( isQuadratic ) {
1213       aResVec[SMDSEntity_Node     ] = 2 * ( layerPositions.size() - 1 ) + 1;
1214       aResVec[SMDSEntity_Quad_Edge] = layerPositions.size() - 1;
1215     }
1216     else {
1217       aResVec[SMDSEntity_Node] = layerPositions.size() - 2;
1218       aResVec[SMDSEntity_Edge] = layerPositions.size() - 1;
1219     }
1220     sm = aMesh.GetSubMesh( linSide1->Edge(0) );
1221     aResMap[sm] = aResVec;
1222     if ( linSide2 )
1223     {
1224       sm = aMesh.GetSubMesh( linSide2->Edge(0) );
1225       aResMap[sm] = aResVec;
1226     }
1227   }
1228
1229   return true;
1230 }
1231
1232 //================================================================================
1233 /*!
1234  * \brief Return true if the algorithm can compute mesh on this shape
1235  */
1236 //================================================================================
1237
1238 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1239 {
1240   int nbFoundFaces = 0;
1241   for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
1242   {
1243     StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1244     int nbSides = analyseFace( exp.Current(), NULL, circSide, linSide1, linSide2, NULL );
1245     bool ok = ( 0 < nbSides && nbSides <= 3 &&
1246                 isCornerInsideCircle( circSide, linSide1, linSide2 ));
1247     if( toCheckAll  && !ok ) return false;
1248     if( !toCheckAll && ok  ) return true;
1249   }
1250   if( toCheckAll && nbFoundFaces != 0 ) return true;
1251   return false;
1252 };