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