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   {
 762     for ( TopExp_Explorer e( faceSubMesh->GetSubShape(), TopAbs_EDGE ); e.More(); e.Next() )
 763     {
 764       TEdgeMarker::markEdge( TopoDS::Edge( e.Current() ), faceSubMesh );
 765     }
 766   }
 767 }
 768
 769 //=======================================================================
 770 //function : Compute
 771 //purpose  :
 772 //=======================================================================
 773
 774 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 775                                                const TopoDS_Shape& aShape)
 776 {
 777   SMESH_MesherHelper helper( aMesh );
 778   StdMeshers_Quadrangle_2D::myHelper     = & helper;
 779   StdMeshers_Quadrangle_2D::myNeedSmooth = false;
 780   StdMeshers_Quadrangle_2D::myCheckOri   = false;
 781   StdMeshers_Quadrangle_2D::myQuadList.clear();
 782
 783   myHelper->SetSubShape( aShape );
 784   myHelper->SetElementsOnShape( true );
 785
 786   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
 787   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
 788   if( nbSides > 3 || nbSides < 1 )
 789     return error("The face must be a full ellipse or a part of ellipse (i.e. the number "
 790                  "of edges is less or equal to 3 and one of them is an ellipse curve)");
 791
 792   // get not yet computed EDGEs
 793   list< TopoDS_Edge > emptyEdges;
 794   for ( TopExp_Explorer e( aShape, TopAbs_EDGE ); e.More(); e.Next() )
 795   {
 796     if ( aMesh.GetSubMesh( e.Current() )->IsEmpty() )
 797       emptyEdges.push_back( TopoDS::Edge( e.Current() ));
 798   }
 799
 800   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
 801
 802   if ( !algo1d->ComputeCircularEdge( aMesh, circSide ))
 803     return error( algo1d->GetComputeError() );
 804
 805
 806   TopoDS_Face F = TopoDS::Face(aShape);
 807   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 808
 809   myHelper->IsQuadraticSubMesh( aShape );
 810
 811   vector< double > layerPositions; // [0,1]
 812
 813   const SMDS_MeshNode* centerNode = 0;
 814   gp_Pnt2d             centerUV(0,0);
 815
 816   // ------------------------------------------------------------------------------------------
 817   if ( nbSides == 1 ) // full ellipse
 818   {
 819     const SMDS_MeshNode* circNode;
 820     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
 821
 822     StdMeshers_FaceSidePtr tmpSide =
 823       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
 824
 825     if ( !computeLayerPositions( tmpSide, layerPositions ))
 826       return false;
 827
 828     UVPtStructVec nodes( layerPositions.size() );
 829     nodes[0].node = circNode;
 830     for ( size_t i = 0; i < layerPositions.size(); ++i )
 831     {
 832       gp_Pnt2d uv = tmpSide->Value2d( layerPositions[i] );
 833       gp_Pnt  xyz = S->Value( uv.X(), uv.Y() );
 834
 835       nodes[ i ].SetUV( uv.XY() );
 836       nodes[ i ].normParam = layerPositions[i];
 837       if ( i )
 838         nodes[ i ].node = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z(), 0, uv.X(), uv.Y() );
 839     }
 840
 841     linSide1 = StdMeshers_FaceSide::New( nodes );
 842     linSide2 = StdMeshers_FaceSide::New( nodes );
 843
 844     centerNode = nodes.back().node;
 845     centerUV   = nodes.back().UV();
 846   }
 847   // ------------------------------------------------------------------------------------------
 848   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
 849   {
 850     // full ellipse with an internal radial side
 851
 852     // eliminate INTERNAL orientation
 853     list< TopoDS_Edge > edges;
 854     for ( int iE = 0; iE < linSide1->NbEdges(); ++iE )
 855     {
 856       edges.push_back( linSide1->Edge( iE ));
 857       edges.back().Orientation( TopAbs_FORWARD );
 858     }
 859
 860     // orient the internal side
 861     bool isVIn0Shared = false;
 862     TopoDS_Vertex vIn0 = myHelper->IthVertex( 0, edges.front() );
 863     for ( int iE = 0; iE < circSide->NbEdges() && !isVIn0Shared; ++iE )
 864       isVIn0Shared = vIn0.IsSame( circSide->FirstVertex( iE ));
 865
 866     linSide1 = StdMeshers_FaceSide::New( F, edges, &aMesh,
 867                                          /*isFrw=*/isVIn0Shared, /*skipMedium=*/true, myHelper );
 868
 869     int nbMeshedEdges;
 870     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges ))
 871       return false;
 872
 873     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 874     UVPtStructVec nodes;
 875     if ( nbMeshedEdges != linSide1->NbEdges() )
 876       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 877     else
 878       nodes = linSide1->GetUVPtStruct();
 879
 880     linSide1 = StdMeshers_FaceSide::New( nodes );
 881     linSide2 = StdMeshers_FaceSide::New( nodes );
 882
 883     centerNode = nodes.back().node;
 884     centerUV   = nodes.back().UV();
 885   }
 886   // ------------------------------------------------------------------------------------------
 887   else if ( nbSides == 2 )
 888   {
 889     // find positions of layers for the first half of linSide1
 890     int nbMeshedEdges;
 891     if ( !computeLayerPositions( linSide1, layerPositions, &nbMeshedEdges, /*useHalf=*/true ))
 892       return false;
 893
 894     // make positions for the whole linSide1
 895     for ( size_t i = 0; i < layerPositions.size(); ++i )
 896     {
 897       layerPositions[i] *= 0.5;
 898     }
 899     layerPositions.reserve( layerPositions.size() * 2 );
 900     for ( int nb = layerPositions.size()-1; nb > 0; --nb )
 901       layerPositions.push_back( layerPositions.back() + layerPositions[nb] - layerPositions[nb-1] );
 902
 903     // merge existing nodes with new nodes at layerPositions into a UVPtStructVec
 904     UVPtStructVec nodes;
 905     if ( nbMeshedEdges != linSide1->NbEdges() )
 906       makeMissingMesh( linSide1, nodes, layerPositions, myHelper );
 907     else
 908       nodes = linSide1->GetUVPtStruct();
 909
 910     // find a central node
 911     size_t i = 0;
 912     while ( nodes[ i ].normParam < 0.5 ) ++i;
 913     if (( nodes[ i ].normParam - 0.5 ) > ( 0.5 - nodes[ i-1 ].normParam )) --i;
 914
 915     // distribute nodes between two linear sides
 916     UVPtStructVec nodes2( nodes.rbegin(), nodes.rbegin() + nodes.size() - i );
 917     nodes.resize( i + 1 );
 918
 919     linSide1 = StdMeshers_FaceSide::New( nodes );
 920     linSide2 = StdMeshers_FaceSide::New( nodes2 );
 921
 922     centerNode = nodes.back().node;
 923     centerUV   = nodes.back().UV();
 924   }
 925   // ------------------------------------------------------------------------------------------
 926   else // nbSides == 3
 927   {
 928     // one curve must be a part of ellipse and 2 other curves must be segments of line
 929
 930     int nbMeshedEdges1, nbMeshedEdges2;
 931     vector< double > layerPositions2;
 932     bool ok1 = computeLayerPositions( linSide1, layerPositions,  &nbMeshedEdges1 );
 933     bool ok2 = computeLayerPositions( linSide2, layerPositions2, &nbMeshedEdges2 );
 934     if ( !ok1 && !ok2 )
 935       return false;
 936
 937     bool linSide1Computed = ( nbMeshedEdges1 == linSide1->NbEdges() );
 938     bool linSide2Computed = ( nbMeshedEdges2 == linSide2->NbEdges() );
 939
 940     UVPtStructVec nodes;
 941
 942     if ( linSide1Computed && !linSide2Computed )
 943     {
 944       // use layer positions of linSide1 to mesh linSide2
 945       makeMissingMesh( linSide2, nodes, layerPositions, myHelper );
 946       linSide2 = StdMeshers_FaceSide::New( nodes );
 947     }
 948     else if ( linSide2Computed && !linSide1Computed )
 949     {
 950       // use layer positions of linSide2 to mesh linSide1
 951       makeMissingMesh( linSide1, nodes, layerPositions2, myHelper );
 952       linSide1 = StdMeshers_FaceSide::New( nodes );
 953     }
 954     else if ( !linSide2Computed && !linSide1Computed )
 955     {
 956       // use layer positions of a longer side to mesh the shorter side
 957       vector< double >& lp =
 958         ( linSide1->Length() > linSide2->Length() ) ? layerPositions : layerPositions2;
 959
 960       makeMissingMesh( linSide1, nodes, lp, myHelper );
 961       linSide1 = StdMeshers_FaceSide::New( nodes );
 962       makeMissingMesh( linSide2, nodes, lp, myHelper );
 963       linSide2 = StdMeshers_FaceSide::New( nodes );
 964     }
 965
 966     const UVPtStructVec& nodes2 = linSide2->GetUVPtStruct();
 967     centerNode = nodes2.back().node;
 968     centerUV   = nodes2.back().UV();
 969   }
 970
 971   list< TopoDS_Edge >::iterator ee = emptyEdges.begin();
 972   for ( ; ee != emptyEdges.end(); ++ee )
 973     TEdgeMarker::markEdge( *ee, aMesh.GetSubMesh( F ));
 974
 975   circSide->GetUVPtStruct(); // let sides take into account just computed nodes
 976   linSide1->GetUVPtStruct();
 977   linSide2->GetUVPtStruct();
 978
 979   FaceQuadStruct::Ptr quad( new FaceQuadStruct );
 980   quad->side.resize( 4 );
 981   quad->side[0] = circSide;
 982   quad->side[1] = linSide1;
 983   quad->side[2] = StdMeshers_FaceSide::New( circSide.get(), centerNode, &centerUV );
 984   quad->side[3] = linSide2;
 985
 986   myQuadList.push_back( quad );
 987
 988   // create quadrangles
 989   bool ok;
 990   if ( linSide1->NbPoints() == linSide2->NbPoints() )
 991     ok = StdMeshers_Quadrangle_2D::computeQuadDominant( aMesh, F, quad );
 992   else
 993     ok = StdMeshers_Quadrangle_2D::computeTriangles( aMesh, F, quad );
 994
 995   StdMeshers_Quadrangle_2D::myHelper = 0;
 996
 997   return ok;
 998 }
 999
1000 //================================================================================
1001 /*!
1002  * \brief Compute nodes on the radial edge
1003  * \retval int - nb of segments on the linSide
1004  */
1005 //================================================================================
1006
1007 int StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(StdMeshers_FaceSidePtr linSide,
1008                                                             vector< double > &     positions,
1009                                                             int*                   nbMeshedEdges,
1010                                                             bool                   useHalf)
1011 {
1012   // First, try to compute positions of layers
1013
1014   positions.clear();
1015
1016   SMESH_Mesh * mesh = myHelper->GetMesh();
1017
1018   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
1019   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
1020
1021   if ( !hyp1D && !nbLayers )
1022   {
1023     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
1024     // find a hyp usable by TNodeDistributor
1025     TopoDS_Shape edge = linSide->Edge(0);
1026     const SMESH_HypoFilter* hypKind =
1027       TNodeDistributor::GetDistributor(*mesh)->GetCompatibleHypoFilter(/*ignoreAux=*/true);
1028     hyp1D = mesh->GetHypothesis( edge, *hypKind, /*fromAncestors=*/true);
1029   }
1030   if ( hyp1D ) // try to compute with hyp1D
1031   {
1032     BRepAdaptor_CompCurve* curve = linSide->GetCurve3d();
1033     SMESHUtils::Deleter< BRepAdaptor_CompCurve > delCurve( curve );
1034     double f = curve->FirstParameter();
1035     double l = curve->LastParameter();
1036
1037     if ( useHalf )
1038       l = 0.5 * ( f + l ); // first part of linSide is used
1039
1040     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( positions, linSide->Edge(0),
1041                                                             *curve, f, l, *mesh, hyp1D ))
1042     {
1043       if ( myDistributionHypo ) { // bad hyp assigned
1044         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1045       }
1046       else {
1047         // bad hyp found, its Ok, lets try with default nb of segments
1048       }
1049     }
1050   }
1051
1052   if ( positions.empty() ) // try to use nb of layers
1053   {
1054     if ( !nbLayers )
1055       nbLayers = _gen->GetDefaultNbSegments();
1056
1057     if ( nbLayers )
1058     {
1059       positions.resize( nbLayers + 1 );
1060       for ( int z = 0; z < nbLayers; ++z )
1061         positions[ z ] = double( z )/ double( nbLayers );
1062       positions.back() = 1;
1063     }
1064   }
1065
1066   // Second, check presence of a mesh built by other algo on linSide
1067
1068   int nbEdgesComputed = 0;
1069   for ( int i = 0; i < linSide->NbEdges(); ++i )
1070   {
1071     nbEdgesComputed += ( !mesh->GetSubMesh( linSide->Edge(i))->IsEmpty() );
1072   }
1073
1074   if ( nbEdgesComputed == linSide->NbEdges() )
1075   {
1076     const UVPtStructVec& points = linSide->GetUVPtStruct();
1077     if ( points.size() >= 2 )
1078     {
1079       positions.resize( points.size() );
1080       for ( size_t i = 0; i < points.size(); ++i )
1081         positions[ i ] = points[i].normParam;
1082     }
1083   }
1084
1085   if ( nbMeshedEdges ) *nbMeshedEdges = nbEdgesComputed;
1086
1087   return positions.size();
1088 }
1089
1090
1091 //=======================================================================
1092 //function : Evaluate
1093 //purpose  :
1094 //=======================================================================
1095
1096 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh&         aMesh,
1097                                                 const TopoDS_Shape& aShape,
1098                                                 MapShapeNbElems&    aResMap)
1099 {
1100   if( aShape.ShapeType() != TopAbs_FACE ) {
1101     return false;
1102   }
1103   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1104   if( aResMap.count(sm) )
1105     return false;
1106
1107   vector<int>& aResVec =
1108     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1109
1110   myHelper = new SMESH_MesherHelper( aMesh );
1111   myHelper->SetSubShape( aShape );
1112   SMESHUtils::Deleter<SMESH_MesherHelper> helperDeleter( myHelper );
1113
1114   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1115
1116   StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1117   int nbSides = analyseFace( aShape, &aMesh, circSide, linSide1, linSide2, myHelper );
1118   if( nbSides > 3 || nbSides < 1 )
1119     return false;
1120
1121   if ( algo1d->EvaluateCircularEdge( aMesh, circSide, aResMap ))
1122     return false;
1123
1124   vector< double > layerPositions; // [0,1]
1125
1126   // ------------------------------------------------------------------------------------------
1127   if ( nbSides == 1 )
1128   {
1129     const TopoDS_Face& F = TopoDS::Face( aShape );
1130
1131     const SMDS_MeshNode* circNode;
1132     TopoDS_Edge linEdge = makeEdgeToCenter( circSide, F, circNode );
1133
1134     StdMeshers_FaceSidePtr tmpSide =
1135       StdMeshers_FaceSide::New( F, linEdge, &aMesh, /*isFrw=*/true, /*skipMedium=*/true, myHelper );
1136
1137     if ( !computeLayerPositions( tmpSide, layerPositions ))
1138       return false;
1139   }
1140   // ------------------------------------------------------------------------------------------
1141   else if ( nbSides == 2 && linSide1->Edge(0).Orientation() == TopAbs_INTERNAL )
1142   {
1143     if ( !computeLayerPositions( linSide1, layerPositions ))
1144       return false;
1145   }
1146   // ------------------------------------------------------------------------------------------
1147   else if ( nbSides == 2 )
1148   {
1149     // find positions of layers for the first half of linSide1
1150     if ( !computeLayerPositions( linSide1, layerPositions, 0, /*useHalf=*/true ))
1151       return false;
1152   }
1153   // ------------------------------------------------------------------------------------------
1154   else // nbSides == 3
1155   {
1156     if ( !computeLayerPositions(( linSide1->Length() > linSide2->Length() ) ? linSide1 : linSide2,
1157                                  layerPositions ))
1158       return false;
1159   }
1160
1161   bool isQuadratic = false;
1162   for ( TopExp_Explorer edge( aShape, TopAbs_EDGE ); edge.More() &&  !isQuadratic ; edge.Next() )
1163   {
1164     sm = aMesh.GetSubMesh( edge.Current() );
1165     vector<int>& nbElems = aResMap[ sm ];
1166     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1167       isQuadratic = nbElems[ SMDSEntity_Quad_Edge ];
1168   }
1169
1170   int nbCircSegments = 0;
1171   for ( int iE = 0; iE < circSide->NbEdges(); ++iE )
1172   {
1173     sm = aMesh.GetSubMesh( circSide->Edge( iE ));
1174     vector<int>& nbElems = aResMap[ sm ];
1175     if ( SMDSEntity_Quad_Edge < (int) nbElems.size() )
1176       nbCircSegments += ( nbElems[ SMDSEntity_Edge ] + nbElems[ SMDSEntity_Quad_Edge ]);
1177   }
1178
1179   int nbQuads = nbCircSegments * ( layerPositions.size() - 1 );
1180   int nbTria  = nbCircSegments;
1181   int nbNodes = ( nbCircSegments - 1 ) * ( layerPositions.size() - 2 );
1182   if ( isQuadratic )
1183   {
1184     nbNodes += (( nbCircSegments - 1 ) * ( layerPositions.size() - 1 ) + // radial
1185                 ( nbCircSegments     ) * ( layerPositions.size() - 2 )); // circular
1186     aResVec[SMDSEntity_Quad_Triangle  ] = nbTria;
1187     aResVec[SMDSEntity_Quad_Quadrangle] = nbQuads;
1188   }
1189   else
1190   {
1191     aResVec[SMDSEntity_Triangle  ] = nbTria;
1192     aResVec[SMDSEntity_Quadrangle] = nbQuads;
1193   }
1194   aResVec[SMDSEntity_Node] = nbNodes;
1195
1196   if ( linSide1 )
1197   {
1198     // evaluation for linSides
1199     vector<int> aResVec(SMDSEntity_Last, 0);
1200     if ( isQuadratic ) {
1201       aResVec[SMDSEntity_Node     ] = 2 * ( layerPositions.size() - 1 ) + 1;
1202       aResVec[SMDSEntity_Quad_Edge] = layerPositions.size() - 1;
1203     }
1204     else {
1205       aResVec[SMDSEntity_Node] = layerPositions.size() - 2;
1206       aResVec[SMDSEntity_Edge] = layerPositions.size() - 1;
1207     }
1208     sm = aMesh.GetSubMesh( linSide1->Edge(0) );
1209     aResMap[sm] = aResVec;
1210     if ( linSide2 )
1211     {
1212       sm = aMesh.GetSubMesh( linSide2->Edge(0) );
1213       aResMap[sm] = aResVec;
1214     }
1215   }
1216
1217   return true;
1218 }
1219
1220 //================================================================================
1221 /*!
1222  * \brief Return true if the algorithm can compute mesh on this shape
1223  */
1224 //================================================================================
1225
1226 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1227 {
1228   int nbFoundFaces = 0;
1229   for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
1230   {
1231     StdMeshers_FaceSidePtr circSide, linSide1, linSide2;
1232     int nbSides = analyseFace( exp.Current(), NULL, circSide, linSide1, linSide2, NULL );
1233     bool ok = ( 0 < nbSides && nbSides <= 3 &&
1234                 isCornerInsideCircle( circSide, linSide1, linSide2 ));
1235     if( toCheckAll  && !ok ) return false;
1236     if( !toCheckAll && ok  ) return true;
1237   }
1238   if( toCheckAll && nbFoundFaces != 0 ) return true;
1239   return false;
1240 };