src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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