src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

   1 // Copyright (C) 2007-2014  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 //  SMESH SMESH : implementaion of SMESH idl descriptions
  21 // File      : StdMeshers_RadialQuadrangle_1D2D.cxx
  22 // Module    : SMESH
  23
  24 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
  25
  26 #include "StdMeshers_NumberOfLayers.hxx"
  27 #include "StdMeshers_LayerDistribution.hxx"
  28 #include "StdMeshers_Regular_1D.hxx"
  29 #include "StdMeshers_NumberOfSegments.hxx"
  30
  31 #include "SMDS_MeshNode.hxx"
  32 #include "SMESHDS_SubMesh.hxx"
  33 #include "SMESH_Gen.hxx"
  34 #include "SMESH_HypoFilter.hxx"
  35 #include "SMESH_Mesh.hxx"
  36 #include "SMESH_MesherHelper.hxx"
  37 #include "SMESH_subMesh.hxx"
  38 #include "SMESH_subMeshEventListener.hxx"
  39
  40 #include "utilities.h"
  41
  42 #include <BRepAdaptor_Curve.hxx>
  43 #include <BRepBuilderAPI_MakeEdge.hxx>
  44 #include <BRep_Tool.hxx>
  45 #include <GeomAPI_ProjectPointOnSurf.hxx>
  46 #include <Geom_Circle.hxx>
  47 #include <Geom_Line.hxx>
  48 #include <Geom_TrimmedCurve.hxx>
  49 #include <TColgp_SequenceOfPnt.hxx>
  50 #include <TColgp_SequenceOfPnt2d.hxx>
  51 #include <TopExp.hxx>
  52 #include <TopExp_Explorer.hxx>
  53 #include <TopTools_ListIteratorOfListOfShape.hxx>
  54 #include <TopoDS.hxx>
  55
  56
  57 using namespace std;
  58
  59 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  60 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  61
  62
  63 //=======================================================================
  64 //function : StdMeshers_RadialQuadrangle_1D2D
  65 //purpose  :
  66 //=======================================================================
  67
  68 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
  69                                                                    int studyId,
  70                                                                    SMESH_Gen* gen)
  71   :SMESH_2D_Algo(hypId, studyId, gen)
  72 {
  73   _name = "RadialQuadrangle_1D2D";
  74   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
  75
  76   _compatibleHypothesis.push_back("LayerDistribution2D");
  77   _compatibleHypothesis.push_back("NumberOfLayers2D");
  78   myNbLayerHypo = 0;
  79   myDistributionHypo = 0;
  80   _requireDiscreteBoundary = false;
  81   _supportSubmeshes = true;
  82 }
  83
  84
  85 //================================================================================
  86 /*!
  87  * \brief Destructor
  88  */
  89 //================================================================================
  90
  91 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
  92 {}
  93
  94
  95 //=======================================================================
  96 //function : CheckHypothesis
  97 //purpose  :
  98 //=======================================================================
  99
 100 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
 101                            (SMESH_Mesh&                          aMesh,
 102                             const TopoDS_Shape&                  aShape,
 103                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
 104 {
 105   // check aShape
 106   myNbLayerHypo = 0;
 107   myDistributionHypo = 0;
 108
 109   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 110
 111   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 112   if ( hyps.size() == 0 ) {
 113     aStatus = SMESH_Hypothesis::HYP_OK;
 114     return true;  // can work with no hypothesis
 115   }
 116
 117   if ( hyps.size() > 1 ) {
 118     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 119     return false;
 120   }
 121
 122   const SMESHDS_Hypothesis *theHyp = hyps.front();
 123
 124   string hypName = theHyp->GetName();
 125
 126   if (hypName == "NumberOfLayers2D") {
 127     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 128     aStatus = SMESH_Hypothesis::HYP_OK;
 129     return true;
 130   }
 131   if (hypName == "LayerDistribution2D") {
 132     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 133     aStatus = SMESH_Hypothesis::HYP_OK;
 134     return true;
 135   }
 136   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 137   return true;
 138 }
 139
 140 namespace
 141 {
 142   // ------------------------------------------------------------------------------
 143   /*!
 144    * \brief Listener used to mark edges meshed by StdMeshers_RadialQuadrangle_1D2D
 145    */
 146   class TEdgeMarker : public SMESH_subMeshEventListener
 147   {
 148     TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
 149                                               "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
 150   public:
 151     //!<  Return static listener
 152     static SMESH_subMeshEventListener* getListener()
 153     {
 154       static TEdgeMarker theEdgeMarker;
 155       return &theEdgeMarker;
 156     }
 157     //! Clear face sumbesh if something happens on edges
 158     void ProcessEvent(const int          event,
 159                       const int          eventType,
 160                       SMESH_subMesh*     edgeSubMesh,
 161                       EventListenerData* data,
 162                       const SMESH_Hypothesis*  /*hyp*/)
 163     {
 164       if ( data && !data->mySubMeshes.empty() && eventType == SMESH_subMesh::ALGO_EVENT)
 165       {
 166         ASSERT( data->mySubMeshes.front() != edgeSubMesh );
 167         SMESH_subMesh* faceSubMesh = data->mySubMeshes.front();
 168         faceSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
 169       }
 170     }
 171   };
 172
 173   // ------------------------------------------------------------------------------
 174   /*!
 175    * \brief Mark an edge as computed by StdMeshers_RadialQuadrangle_1D2D
 176    */
 177   void markEdgeAsComputedByMe(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
 178   {
 179     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
 180     {
 181       if ( !edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
 182         faceSubMesh->SetEventListener( TEdgeMarker::getListener(),
 183                                        SMESH_subMeshEventListenerData::MakeData(faceSubMesh),
 184                                        edgeSM);
 185     }
 186   }
 187   // ------------------------------------------------------------------------------
 188   /*!
 189    * \brief Return true if a radial edge was meshed with StdMeshers_RadialQuadrangle_1D2D with
 190    * the same radial distribution
 191    */
 192 //   bool isEdgeCompatiballyMeshed(const TopoDS_Edge& edge, SMESH_subMesh* faceSubMesh)
 193 //   {
 194 //     if ( SMESH_subMesh* edgeSM = faceSubMesh->GetFather()->GetSubMeshContaining( edge ))
 195 //     {
 196 //       if ( SMESH_subMeshEventListenerData* otherFaceData =
 197 //            edgeSM->GetEventListenerData( TEdgeMarker::getListener() ))
 198 //       {
 199 //         // compare hypothesis aplied to two disk faces sharing radial edges
 200 //         SMESH_Mesh& mesh = *faceSubMesh->GetFather();
 201 //         SMESH_Algo* radialQuadAlgo = mesh.GetGen()->GetAlgo(mesh, faceSubMesh->GetSubShape() );
 202 //         SMESH_subMesh* otherFaceSubMesh = otherFaceData->mySubMeshes.front();
 203 //         list <const SMESHDS_Hypothesis *> hyps1 =
 204 //           radialQuadAlgo->GetUsedHypothesis( mesh, faceSubMesh->GetSubShape());
 205 //         list <const SMESHDS_Hypothesis *> hyps2 =
 206 //           radialQuadAlgo->GetUsedHypothesis( mesh, otherFaceSubMesh->GetSubShape());
 207 //         if( hyps1.empty() && hyps2.empty() )
 208 //           return true; // defaul hyps
 209 //         if ( hyps1.size() != hyps2.size() )
 210 //           return false;
 211 //         return *hyps1.front() == *hyps2.front();
 212 //       }
 213 //     }
 214 //     return false;
 215 //   }
 216
 217   //================================================================================
 218   /*!
 219    * \brief Return base curve of the edge and extremum parameters
 220    */
 221   //================================================================================
 222
 223   Handle(Geom_Curve) getCurve(const TopoDS_Edge& edge, double* f=0, double* l=0)
 224   {
 225     Handle(Geom_Curve) C;
 226     if ( !edge.IsNull() )
 227     {
 228       double first = 0., last = 0.;
 229       C = BRep_Tool::Curve(edge, first, last);
 230       if ( !C.IsNull() )
 231       {
 232         Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C);
 233         while( !tc.IsNull() ) {
 234           C = tc->BasisCurve();
 235           tc = Handle(Geom_TrimmedCurve)::DownCast(C);
 236         }
 237         if ( f ) *f = first;
 238         if ( l ) *l = last;
 239       }
 240     }
 241     return C;
 242   }
 243
 244   //================================================================================
 245   /*!
 246    * \brief Return edges of the face
 247    *  \retval int - nb of edges
 248    */
 249   //================================================================================
 250
 251   int analyseFace(const TopoDS_Shape& face,
 252                   TopoDS_Edge&        CircEdge,
 253                   TopoDS_Edge&        LinEdge1,
 254                   TopoDS_Edge&        LinEdge2)
 255   {
 256     CircEdge.Nullify(); LinEdge1.Nullify(); LinEdge2.Nullify();
 257     int nbe = 0;
 258
 259     for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next(), ++nbe )
 260     {
 261       const TopoDS_Edge& E = TopoDS::Edge( exp.Current() );
 262       double f,l;
 263       Handle(Geom_Curve) C = getCurve(E,&f,&l);
 264       if ( !C.IsNull() )
 265       {
 266         if ( C->IsKind( STANDARD_TYPE(Geom_Circle)))
 267         {
 268           if ( CircEdge.IsNull() )
 269             CircEdge = E;
 270           else
 271             return 0;
 272         }
 273         else if ( LinEdge1.IsNull() )
 274           LinEdge1 = E;
 275         else
 276           LinEdge2 = E;
 277       }
 278     }
 279     return nbe;
 280   }
 281
 282 //================================================================================
 283 //================================================================================
 284 /*!
 285  * \brief Class computing layers distribution using data of
 286  *        StdMeshers_LayerDistribution hypothesis
 287  */
 288 //================================================================================
 289 //================================================================================
 290
 291 class TNodeDistributor: public StdMeshers_Regular_1D
 292 {
 293   list <const SMESHDS_Hypothesis *> myUsedHyps;
 294 public:
 295   // -----------------------------------------------------------------------------
 296   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 297   {
 298     const int myID = -1001;
 299     TNodeDistributor* myHyp = dynamic_cast<TNodeDistributor*>( aMesh.GetHypothesis( myID ));
 300     if ( !myHyp )
 301       myHyp = new TNodeDistributor( myID, 0, aMesh.GetGen() );
 302     return myHyp;
 303   }
 304   // -----------------------------------------------------------------------------
 305   //! Computes distribution of nodes on a straight line ending at pIn and pOut
 306   bool Compute( vector< double > &      positions,
 307                 gp_Pnt                  pIn,
 308                 gp_Pnt                  pOut,
 309                 SMESH_Mesh&             aMesh,
 310                 const SMESH_Hypothesis* hyp1d)
 311   {
 312     if ( !hyp1d ) return error( "Invalid LayerDistribution hypothesis");
 313
 314     double len = pIn.Distance( pOut );
 315     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
 316
 317     myUsedHyps.clear();
 318     myUsedHyps.push_back( hyp1d );
 319
 320     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
 321     SMESH_Hypothesis::Hypothesis_Status aStatus;
 322     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 323       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 324                     "with LayerDistribution hypothesis");
 325
 326     BRepAdaptor_Curve C3D(edge);
 327     double f = C3D.FirstParameter(), l = C3D.LastParameter();
 328     list< double > params;
 329     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
 330       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 331
 332     positions.clear();
 333     positions.reserve( params.size() );
 334     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 335       positions.push_back( *itU / len );
 336     return true;
 337   }
 338   // -----------------------------------------------------------------------------
 339   //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
 340   bool ComputeCircularEdge(SMESH_Mesh&         aMesh,
 341                            const TopoDS_Edge& anEdge)
 342   {
 343     _gen->Compute( aMesh, anEdge);
 344     SMESH_subMesh *sm = aMesh.GetSubMesh(anEdge);
 345     if ( sm->GetComputeState() != SMESH_subMesh::COMPUTE_OK)
 346     {
 347       // find any 1d hyp assigned (there can be a hyp w/o algo)
 348       myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
 349       Hypothesis_Status aStatus;
 350       if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
 351       {
 352         // no valid 1d hyp assigned, use default nb of segments
 353         _hypType                    = NB_SEGMENTS;
 354         _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 355         _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 356       }
 357       return StdMeshers_Regular_1D::Compute( aMesh, anEdge );
 358     }
 359     return true;
 360   }
 361   // -----------------------------------------------------------------------------
 362   //! Make mesh on an adge using assigned 1d hyp or defaut nb of segments
 363   bool EvaluateCircularEdge(SMESH_Mesh&        aMesh,
 364                             const TopoDS_Edge& anEdge,
 365                             MapShapeNbElems&   aResMap)
 366   {
 367     _gen->Evaluate( aMesh, anEdge, aResMap );
 368     if ( aResMap.count( aMesh.GetSubMesh( anEdge )))
 369       return true;
 370
 371     // find any 1d hyp assigned
 372     myUsedHyps = SMESH_Algo::GetUsedHypothesis(aMesh, anEdge, /*ignoreAux=*/true);
 373     Hypothesis_Status aStatus;
 374     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, anEdge, aStatus ))
 375     {
 376       // no valid 1d hyp assigned, use default nb of segments
 377       _hypType                    = NB_SEGMENTS;
 378       _ivalue[ DISTR_TYPE_IND ]   = StdMeshers_NumberOfSegments::DT_Regular;
 379       _ivalue[ NB_SEGMENTS_IND  ] = _gen->GetDefaultNbSegments();
 380     }
 381     return StdMeshers_Regular_1D::Evaluate( aMesh, anEdge, aResMap );
 382   }
 383 protected:
 384   // -----------------------------------------------------------------------------
 385   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
 386     : StdMeshers_Regular_1D( hypId, studyId, gen)
 387   {
 388   }
 389   // -----------------------------------------------------------------------------
 390   virtual const list <const SMESHDS_Hypothesis *> &
 391     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 392   {
 393     return myUsedHyps;
 394   }
 395   // -----------------------------------------------------------------------------
 396 };
 397 }
 398
 399 //=======================================================================
 400 /*!
 401  * \brief Allow algo to do something after persistent restoration
 402  * \param subMesh - restored submesh
 403  *
 404  * call markEdgeAsComputedByMe()
 405  */
 406 //=======================================================================
 407
 408 void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMesh)
 409 {
 410   if ( !faceSubMesh->IsEmpty() )
 411   {
 412     TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 413     analyseFace( faceSubMesh->GetSubShape(), CircEdge, LinEdge1, LinEdge2 );
 414     if ( !CircEdge.IsNull() ) markEdgeAsComputedByMe( CircEdge, faceSubMesh );
 415     if ( !LinEdge1.IsNull() ) markEdgeAsComputedByMe( LinEdge1, faceSubMesh );
 416     if ( !LinEdge2.IsNull() ) markEdgeAsComputedByMe( LinEdge2, faceSubMesh );
 417   }
 418 }
 419
 420 //=======================================================================
 421 //function : Compute
 422 //purpose  :
 423 //=======================================================================
 424
 425 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 426                                                const TopoDS_Shape& aShape)
 427 {
 428   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 429
 430   myHelper = new SMESH_MesherHelper( aMesh );
 431   myHelper->IsQuadraticSubMesh( aShape );
 432   // to delete helper at exit from Compute()
 433   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
 434
 435   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
 436
 437   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 438   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
 439   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
 440   if( nbe>3 || nbe < 1 || aCirc.IsNull() )
 441     return error("The face must be a full circle or a part of circle (i.e. the number of edges is less or equal to 3 and one of them is a circle curve)");
 442
 443   gp_Pnt P0, P1;
 444   // points for rotation
 445   TColgp_SequenceOfPnt Points;
 446   // angles for rotation
 447   TColStd_SequenceOfReal Angles;
 448   // Nodes1 and Nodes2 - nodes along radiuses
 449   // CNodes - nodes on circle edge
 450   vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
 451   SMDS_MeshNode * NC;
 452   // parameters edge nodes on face
 453   TColgp_SequenceOfPnt2d Pnts2d1;
 454   gp_Pnt2d PC;
 455
 456   int faceID = meshDS->ShapeToIndex(aShape);
 457   TopoDS_Face F = TopoDS::Face(aShape);
 458   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 459
 460
 461   if(nbe==1)
 462   {
 463     if (!algo1d->ComputeCircularEdge( aMesh, CircEdge ))
 464       return error( algo1d->GetComputeError() );
 465     map< double, const SMDS_MeshNode* > theNodes;
 466     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
 467       return error("Circular edge is incorrectly meshed");
 468
 469     CNodes.clear();
 470     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 471     const SMDS_MeshNode* NF = (*itn).second;
 472     CNodes.push_back( (*itn).second );
 473     double fang = (*itn).first;
 474     if ( itn != theNodes.end() ) {
 475       itn++;
 476       for(; itn != theNodes.end(); itn++ ) {
 477         CNodes.push_back( (*itn).second );
 478         double ang = (*itn).first - fang;
 479         if( ang>M_PI ) ang = ang - 2.*M_PI;
 480         if( ang<-M_PI ) ang = ang + 2.*M_PI;
 481         Angles.Append( ang );
 482       }
 483     }
 484     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 485     P0 = aCirc->Location();
 486
 487     if ( !computeLayerPositions(P0,P1))
 488       return false;
 489
 490     TopoDS_Vertex V1 = myHelper->IthVertex(0, CircEdge );
 491     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
 492
 493     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 494     GeomAPI_ProjectPointOnSurf PPS(P0,S);
 495     double U0,V0;
 496     PPS.Parameters(1,U0,V0);
 497     meshDS->SetNodeOnFace(NC, faceID, U0, V0);
 498     PC = gp_Pnt2d(U0,V0);
 499
 500     gp_Vec aVec(P0,P1);
 501     gp_Vec2d aVec2d(PC,p2dV);
 502     Nodes1.resize( myLayerPositions.size()+1 );
 503     Nodes2.resize( myLayerPositions.size()+1 );
 504     int i = 0;
 505     for(; i<myLayerPositions.size(); i++) {
 506       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 507                 P0.Y() + aVec.Y()*myLayerPositions[i],
 508                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 509       Points.Append(P);
 510       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 511       Nodes1[i] = node;
 512       Nodes2[i] = node;
 513       double U = PC.X() + aVec2d.X()*myLayerPositions[i];
 514       double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
 515       meshDS->SetNodeOnFace( node, faceID, U, V );
 516       Pnts2d1.Append(gp_Pnt2d(U,V));
 517     }
 518     Nodes1[Nodes1.size()-1] = NF;
 519     Nodes2[Nodes1.size()-1] = NF;
 520   }
 521   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL )
 522   {
 523     // one curve must be a half of circle and other curve must be
 524     // a segment of line
 525     double fp, lp;
 526     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
 527     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
 528       // not half of circle
 529       return error(COMPERR_BAD_SHAPE);
 530     }
 531     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
 532     if( aLine.IsNull() ) {
 533       // other curve not line
 534       return error(COMPERR_BAD_SHAPE);
 535     }
 536
 537     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
 538       return error( algo1d->GetComputeError() );
 539     map< double, const SMDS_MeshNode* > theNodes;
 540     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) )
 541       return error("Circular edge is incorrectly meshed");
 542
 543     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 544     CNodes.clear();
 545     CNodes.push_back( itn->second );
 546     double fang = (*itn).first;
 547     itn++;
 548     for(; itn != theNodes.end(); itn++ ) {
 549       CNodes.push_back( (*itn).second );
 550       double ang = (*itn).first - fang;
 551       if( ang>M_PI ) ang = ang - 2.*M_PI;
 552       if( ang<-M_PI ) ang = ang + 2.*M_PI;
 553       Angles.Append( ang );
 554     }
 555     const SMDS_MeshNode* NF = theNodes.begin()->second;
 556     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
 557     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 558     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 559     P0 = aCirc->Location();
 560
 561     bool linEdgeComputed;
 562     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdgeComputed))
 563       return false;
 564
 565     if ( linEdgeComputed )
 566     {
 567       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
 568         return error("Invalid mesh on a straight edge");
 569
 570       Nodes1.resize( myLayerPositions.size()+1 );
 571       Nodes2.resize( myLayerPositions.size()+1 );
 572       vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
 573       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
 574       if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
 575
 576       map< double, const SMDS_MeshNode* >::reverse_iterator ritn = theNodes.rbegin();
 577       itn = theNodes.begin();
 578       for ( int i = Nodes1.size()-1; i > -1; ++itn, ++ritn, --i )
 579       {
 580         (*pNodes1)[i] = ritn->second;
 581         (*pNodes2)[i] =  itn->second;
 582         Points.Prepend( gpXYZ( Nodes1[i]));
 583         Pnts2d1.Prepend( myHelper->GetNodeUV( F, Nodes1[i]));
 584       }
 585       NC = const_cast<SMDS_MeshNode*>( itn->second );
 586       Points.Remove( Nodes1.size() );
 587     }
 588     else
 589     {
 590       gp_Vec aVec(P0,P1);
 591       int edgeID = meshDS->ShapeToIndex(LinEdge1);
 592       // check orientation
 593       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 594       gp_Pnt Ptmp;
 595       Crv->D0(fp,Ptmp);
 596       bool ori = true;
 597       if( P1.Distance(Ptmp) > Precision::Confusion() )
 598         ori = false;
 599       // get UV points for edge
 600       gp_Pnt2d PF,PL;
 601       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 602       PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
 603       gp_Vec2d V2d;
 604       if(ori) V2d = gp_Vec2d(PC,PF);
 605       else V2d = gp_Vec2d(PC,PL);
 606       // add nodes on edge
 607       double cp = (fp+lp)/2;
 608       double dp2 = (lp-fp)/2;
 609       NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 610       meshDS->SetNodeOnEdge(NC, edgeID, cp);
 611       Nodes1.resize( myLayerPositions.size()+1 );
 612       Nodes2.resize( myLayerPositions.size()+1 );
 613       int i = 0;
 614       for(; i<myLayerPositions.size(); i++) {
 615         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 616                   P0.Y() + aVec.Y()*myLayerPositions[i],
 617                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 618         Points.Append(P);
 619         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 620         Nodes1[i] = node;
 621         double param;
 622         if(ori)
 623           param = fp + dp2*(1-myLayerPositions[i]);
 624         else
 625           param = cp + dp2*myLayerPositions[i];
 626         meshDS->SetNodeOnEdge(node, edgeID, param);
 627         P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
 628                     P0.Y() - aVec.Y()*myLayerPositions[i],
 629                     P0.Z() - aVec.Z()*myLayerPositions[i] );
 630         node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 631         Nodes2[i] = node;
 632         if(!ori)
 633           param = fp + dp2*(1-myLayerPositions[i]);
 634         else
 635           param = cp + dp2*myLayerPositions[i];
 636         meshDS->SetNodeOnEdge(node, edgeID, param);
 637         // parameters on face
 638         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 639                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 640         Pnts2d1.Append(P2d);
 641       }
 642       Nodes1[ myLayerPositions.size() ] = NF;
 643       Nodes2[ myLayerPositions.size() ] = NL;
 644       // create 1D elements on edge
 645       vector< const SMDS_MeshNode* > tmpNodes;
 646       tmpNodes.resize(2*Nodes1.size()+1);
 647       for(i=0; i<Nodes2.size(); i++)
 648         tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
 649       tmpNodes[Nodes2.size()] = NC;
 650       for(i=0; i<Nodes1.size(); i++)
 651         tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
 652       for(i=1; i<tmpNodes.size(); i++) {
 653         SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
 654         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 655       }
 656       markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
 657     }
 658   }
 659   else // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
 660   {
 661     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
 662       LinEdge2 = LinEdge1;
 663
 664     // one curve must be a part of circle and other curves must be
 665     // segments of line
 666     double fp, lp;
 667     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
 668     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
 669     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
 670     if( aCirc.IsNull() || aLine1.IsNull() || aLine2.IsNull() )
 671       return error(COMPERR_BAD_SHAPE);
 672
 673     if ( !algo1d->ComputeCircularEdge( aMesh, CircEdge ))
 674       return error( algo1d->GetComputeError() );
 675     map< double, const SMDS_MeshNode* > theNodes;
 676     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes))
 677       return error("Circular edge is incorrectly meshed");
 678
 679     const SMDS_MeshNode* NF = theNodes.begin()->second;
 680     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
 681     CNodes.clear();
 682     CNodes.push_back( NF );
 683     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 684     double fang = (*itn).first;
 685     itn++;
 686     for(; itn != theNodes.end(); itn++ ) {
 687       CNodes.push_back( (*itn).second );
 688       double ang = (*itn).first - fang;
 689       if( ang>M_PI ) ang = ang - 2.*M_PI;
 690       if( ang<-M_PI ) ang = ang + 2.*M_PI;
 691       Angles.Append( ang );
 692     }
 693     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 694     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 695     P0 = aCirc->Location();
 696
 697     // make P1 belong to LinEdge1
 698     TopoDS_Vertex V1 = myHelper->IthVertex( 0, LinEdge1 );
 699     TopoDS_Vertex V2 = myHelper->IthVertex( 1, LinEdge1 );
 700     gp_Pnt PE1 = BRep_Tool::Pnt(V1);
 701     gp_Pnt PE2 = BRep_Tool::Pnt(V2);
 702     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 703         ( P1.Distance(PE2) > Precision::Confusion() ) )
 704       std::swap( LinEdge1, LinEdge2 );
 705
 706     bool linEdge1Computed, linEdge2Computed;
 707     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed))
 708       return false;
 709
 710     Nodes1.resize( myLayerPositions.size()+1 );
 711     Nodes2.resize( myLayerPositions.size()+1 );
 712
 713     // check that both linear edges have same hypotheses
 714     if ( !computeLayerPositions(P0,P2,LinEdge2, &linEdge2Computed))
 715          return false;
 716     if ( Nodes1.size() != myLayerPositions.size()+1 )
 717       return error("Different hypotheses apply to radial edges");
 718
 719     // find the central vertex
 720     TopoDS_Vertex VC = V2;
 721     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 722         ( P2.Distance(PE1) > Precision::Confusion() ) )
 723       VC = V1;
 724     int vertID = meshDS->ShapeToIndex(VC);
 725
 726     // LinEdge1
 727     if ( linEdge1Computed )
 728     {
 729       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
 730         return error("Invalid mesh on a straight edge");
 731
 732       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
 733       NC = const_cast<SMDS_MeshNode*>
 734         ( nodesFromP0ToP1 ? theNodes.begin()->second : theNodes.rbegin()->second );
 735       int i = 0, ir = Nodes1.size()-1;
 736       int * pi = nodesFromP0ToP1 ? &i : &ir;
 737       itn = theNodes.begin();
 738       if ( nodesFromP0ToP1 ) ++itn;
 739       for ( ; i < Nodes1.size(); ++i, --ir, ++itn )
 740       {
 741         Nodes1[*pi] = itn->second;
 742       }
 743       for ( i = 0; i < Nodes1.size()-1; ++i )
 744       {
 745         Points.Append( gpXYZ( Nodes1[i]));
 746         Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
 747       }
 748     }
 749     else
 750     {
 751       int edgeID = meshDS->ShapeToIndex(LinEdge1);
 752       gp_Vec aVec(P0,P1);
 753       // check orientation
 754       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 755       gp_Pnt Ptmp = Crv->Value(fp);
 756       bool ori = false;
 757       if( P1.Distance(Ptmp) > Precision::Confusion() )
 758         ori = true;
 759       // get UV points for edge
 760       gp_Pnt2d PF,PL;
 761       BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 762       gp_Vec2d V2d;
 763       if(ori) {
 764         V2d = gp_Vec2d(PF,PL);
 765         PC = PF;
 766       }
 767       else {
 768         V2d = gp_Vec2d(PL,PF);
 769         PC = PL;
 770       }
 771       NC = const_cast<SMDS_MeshNode*>( VertexNode( VC, meshDS ));
 772       if ( !NC )
 773       {
 774         NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 775         meshDS->SetNodeOnVertex(NC, vertID);
 776       }
 777       double dp = lp-fp;
 778       int i = 0;
 779       for(; i<myLayerPositions.size(); i++) {
 780         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 781                   P0.Y() + aVec.Y()*myLayerPositions[i],
 782                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 783         Points.Append(P);
 784         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 785         Nodes1[i] = node;
 786         double param;
 787         if(!ori)
 788           param = fp + dp*(1-myLayerPositions[i]);
 789         else
 790           param = fp + dp*myLayerPositions[i];
 791         meshDS->SetNodeOnEdge(node, edgeID, param);
 792         // parameters on face
 793         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 794                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 795         Pnts2d1.Append(P2d);
 796       }
 797       Nodes1[ myLayerPositions.size() ] = NF;
 798       // create 1D elements on edge
 799       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
 800       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 801       for(i=1; i<Nodes1.size(); i++) {
 802         ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
 803         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 804       }
 805       if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
 806         Nodes2 = Nodes1;
 807     }
 808     markEdgeAsComputedByMe( LinEdge1, aMesh.GetSubMesh( F ));
 809
 810     // LinEdge2
 811     if ( linEdge2Computed )
 812     {
 813       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge2,true,theNodes))
 814         return error("Invalid mesh on a straight edge");
 815
 816       bool nodesFromP0ToP2 = ( theNodes.rbegin()->second == NL );
 817       int i = 0, ir = Nodes1.size()-1;
 818       int * pi = nodesFromP0ToP2 ? &i : &ir;
 819       itn = theNodes.begin();
 820       if ( nodesFromP0ToP2 ) ++itn;
 821       for ( ; i < Nodes2.size(); ++i, --ir, ++itn )
 822         Nodes2[*pi] = itn->second;
 823     }
 824     else
 825     {
 826       int edgeID = meshDS->ShapeToIndex(LinEdge2);
 827       gp_Vec aVec = gp_Vec(P0,P2);
 828       // check orientation
 829       Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
 830       gp_Pnt Ptmp = Crv->Value(fp);
 831       bool ori = false;
 832       if( P2.Distance(Ptmp) > Precision::Confusion() )
 833         ori = true;
 834       // get UV points for edge
 835       gp_Pnt2d PF,PL;
 836       BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
 837       gp_Vec2d V2d;
 838       if(ori) {
 839         V2d = gp_Vec2d(PF,PL);
 840         PC = PF;
 841       }
 842       else {
 843         V2d = gp_Vec2d(PL,PF);
 844         PC = PL;
 845       }
 846       double dp = lp-fp;
 847       for(int i=0; i<myLayerPositions.size(); i++) {
 848         gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 849                   P0.Y() + aVec.Y()*myLayerPositions[i],
 850                   P0.Z() + aVec.Z()*myLayerPositions[i] );
 851         SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 852         Nodes2[i] = node;
 853         double param;
 854         if(!ori)
 855           param = fp + dp*(1-myLayerPositions[i]);
 856         else
 857           param = fp + dp*myLayerPositions[i];
 858         meshDS->SetNodeOnEdge(node, edgeID, param);
 859         // parameters on face
 860         gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 861                       PC.Y() + V2d.Y()*myLayerPositions[i] );
 862       }
 863       Nodes2[ myLayerPositions.size() ] = NL;
 864       // create 1D elements on edge
 865       SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes2[0] );
 866       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 867       for(int i=1; i<Nodes2.size(); i++) {
 868         ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
 869         if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 870       }
 871     }
 872     markEdgeAsComputedByMe( LinEdge2, aMesh.GetSubMesh( F ));
 873   }
 874   markEdgeAsComputedByMe( CircEdge, aMesh.GetSubMesh( F ));
 875
 876   // orientation
 877   bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
 878   const double angleSign = ( F.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0 );
 879
 880   // create nodes and mesh elements on face
 881   // find axis of rotation
 882   gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
 883   gp_Vec Vec1(P0,P1);
 884   gp_Vec Vec2(P0,P2);
 885   gp_Vec Axis = Vec1.Crossed(Vec2);
 886   // create elements
 887   int i = 1;
 888   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
 889   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
 890   for(; i<Angles.Length(); i++) {
 891     vector< const SMDS_MeshNode* > tmpNodes;
 892     tmpNodes.reserve(Nodes1.size());
 893     gp_Trsf aTrsf;
 894     gp_Ax1 theAxis(P0,gp_Dir(Axis));
 895     aTrsf.SetRotation( theAxis, Angles.Value(i) );
 896     gp_Trsf2d aTrsf2d;
 897     aTrsf2d.SetRotation( PC, Angles.Value(i) * angleSign );
 898     // create nodes
 899     int j = 1;
 900     for(; j<=Points.Length(); j++) {
 901       double cx,cy,cz;
 902       Points.Value(j).Coord( cx, cy, cz );
 903       aTrsf.Transforms( cx, cy, cz );
 904       SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
 905       // find parameters on face
 906       Pnts2d1.Value(j).Coord( cx, cy );
 907       aTrsf2d.Transforms( cx, cy );
 908       // set node on face
 909       meshDS->SetNodeOnFace( node, faceID, cx, cy );
 910       tmpNodes[j-1] = node;
 911     }
 912     // create faces
 913     tmpNodes[Points.Length()] = CNodes[i];
 914     // quad
 915     for(j=0; j<Nodes1.size()-1; j++) {
 916       SMDS_MeshFace* MF;
 917       if(IsForward)
 918         MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
 919                                 Nodes1[j+1], tmpNodes[j+1] );
 920       else
 921         MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
 922                                 Nodes1[j+1], Nodes1[j] );
 923       if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 924     }
 925     // tria
 926     SMDS_MeshFace* MF;
 927     if(IsForward)
 928       MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
 929     else
 930       MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
 931     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 932     for(j=0; j<Nodes1.size(); j++) {
 933       Nodes1[j] = tmpNodes[j];
 934     }
 935   }
 936   // create last faces
 937   // quad
 938   for(i=0; i<Nodes1.size()-1; i++) {
 939     SMDS_MeshFace* MF;
 940     if(IsForward)
 941       MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
 942                               Nodes1[i+1], Nodes2[i+1] );
 943     else
 944       MF = myHelper->AddFace( Nodes2[i],  Nodes2[i+1],
 945                               Nodes1[i+1], Nodes1[i] );
 946     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 947   }
 948   // tria
 949   SMDS_MeshFace* MF;
 950   if(IsForward)
 951     MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
 952   else
 953     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
 954   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 955
 956   return true;
 957 }
 958
 959 //================================================================================
 960 /*!
 961  * \brief Compute positions of nodes on the radial edge
 962   * \retval bool - is a success
 963  */
 964 //================================================================================
 965
 966 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&      p1,
 967                                                              const gp_Pnt&      p2,
 968                                                              const TopoDS_Edge& linEdge,
 969                                                              bool*              linEdgeComputed)
 970 {
 971   // First, try to compute positions of layers
 972
 973   myLayerPositions.clear();
 974
 975   SMESH_Mesh * mesh = myHelper->GetMesh();
 976
 977   const SMESH_Hypothesis* hyp1D = myDistributionHypo ? myDistributionHypo->GetLayerDistribution() : 0;
 978   int                  nbLayers = myNbLayerHypo ? myNbLayerHypo->GetNumberOfLayers() : 0;
 979
 980   if ( !hyp1D && !nbLayers )
 981   {
 982     // No own algo hypotheses assigned, so first try to find any 1D hypothesis.
 983     // We need some edge
 984     TopoDS_Shape edge = linEdge;
 985     if ( edge.IsNull() && !myHelper->GetSubShape().IsNull())
 986       for ( TopExp_Explorer e(myHelper->GetSubShape(), TopAbs_EDGE); e.More(); e.Next())
 987         edge = e.Current();
 988     if ( !edge.IsNull() )
 989     {
 990       // find a hyp usable by TNodeDistributor
 991       SMESH_HypoFilter hypKind;
 992       TNodeDistributor::GetDistributor(*mesh)->InitCompatibleHypoFilter(hypKind,/*ignoreAux=*/1);
 993       hyp1D = mesh->GetHypothesis( edge, hypKind, /*fromAncestors=*/true);
 994     }
 995   }
 996   if ( hyp1D ) // try to compute with hyp1D
 997   {
 998     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) {
 999       if ( myDistributionHypo ) { // bad hyp assigned
1000         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
1001       }
1002       else {
1003         // bad hyp found, its Ok, lets try with default nb of segnents
1004       }
1005     }
1006   }
1007
1008   if ( myLayerPositions.empty() ) // try to use nb of layers
1009   {
1010     if ( !nbLayers )
1011       nbLayers = _gen->GetDefaultNbSegments();
1012
1013     if ( nbLayers )
1014     {
1015       myLayerPositions.resize( nbLayers - 1 );
1016       for ( int z = 1; z < nbLayers; ++z )
1017         myLayerPositions[ z - 1 ] = double( z )/ double( nbLayers );
1018     }
1019   }
1020
1021   // Second, check presence of a mesh built by other algo on linEdge
1022   // and mesh conformity to my hypothesis
1023
1024   bool meshComputed = (!linEdge.IsNull() && !mesh->GetSubMesh(linEdge)->IsEmpty() );
1025   if ( linEdgeComputed ) *linEdgeComputed = meshComputed;
1026
1027   if ( meshComputed )
1028   {
1029     vector< double > nodeParams;
1030     GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams );
1031
1032     // nb of present nodes must be different in cases of 1 and 2 straight edges
1033
1034     TopoDS_Vertex VV[2];
1035     TopExp::Vertices( linEdge, VV[0], VV[1]);
1036     const gp_Pnt* points[] = { &p1, &p2 };
1037     gp_Pnt       vPoints[] = { BRep_Tool::Pnt(VV[0]), BRep_Tool::Pnt(VV[1]) };
1038     const double     tol[] = { BRep_Tool::Tolerance(VV[0]), BRep_Tool::Tolerance(VV[1]) };
1039     bool pointsAreOnVertices = true;
1040     for ( int iP = 0; iP < 2 && pointsAreOnVertices; ++iP )
1041       pointsAreOnVertices = ( points[iP]->Distance( vPoints[0] ) < tol[0] ||
1042                               points[iP]->Distance( vPoints[1] ) < tol[1] );
1043
1044     int nbNodes = nodeParams.size() - 2; // 2 straight edges
1045     if ( !pointsAreOnVertices )
1046       nbNodes = ( nodeParams.size() - 3 ) / 2; // 1 straight edge
1047
1048     if ( myLayerPositions.empty() )
1049     {
1050       myLayerPositions.resize( nbNodes );
1051     }
1052     else if ( myDistributionHypo || myNbLayerHypo )
1053     {
1054       // linEdge is computed by other algo. Check if there is a meshed face
1055       // using nodes on linEdge
1056       bool nodesAreUsed = false;
1057       TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( linEdge );
1058       for ( ; ancestIt.More() && !nodesAreUsed; ancestIt.Next() )
1059         if ( ancestIt.Value().ShapeType() == TopAbs_FACE )
1060           nodesAreUsed = (!mesh->GetSubMesh( ancestIt.Value() )->IsEmpty());
1061       if ( !nodesAreUsed ) {
1062         // rebuild them
1063         mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN );
1064         if ( linEdgeComputed ) *linEdgeComputed = false;
1065       }
1066       else {
1067
1068         if ( myLayerPositions.size() != nbNodes )
1069           return error("Radial edge is meshed by other algorithm");
1070       }
1071     }
1072   }
1073
1074   return !myLayerPositions.empty();
1075 }
1076
1077
1078 //=======================================================================
1079 //function : Evaluate
1080 //purpose  :
1081 //=======================================================================
1082
1083 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
1084                                                 const TopoDS_Shape& aShape,
1085                                                 MapShapeNbElems& aResMap)
1086 {
1087   if( aShape.ShapeType() != TopAbs_FACE ) {
1088     return false;
1089   }
1090   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1091   if( aResMap.count(sm) )
1092     return false;
1093
1094   vector<int>& aResVec =
1095     aResMap.insert( make_pair(sm, vector<int>(SMDSEntity_Last,0))).first->second;
1096
1097   myHelper = new SMESH_MesherHelper( aMesh );
1098   myHelper->SetSubShape( aShape );
1099   auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
1100
1101   TNodeDistributor* algo1d = TNodeDistributor::GetDistributor(aMesh);
1102
1103   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1104   int nbe = analyseFace( aShape, CircEdge, LinEdge1, LinEdge2 );
1105   if( nbe>3 || nbe < 1 || CircEdge.IsNull() )
1106     return false;
1107
1108   Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
1109   if( aCirc.IsNull() )
1110     return error(COMPERR_BAD_SHAPE);
1111
1112   gp_Pnt P0 = aCirc->Location();
1113   gp_Pnt P1 = aCirc->Value(0.);
1114   computeLayerPositions( P0, P1, LinEdge1 );
1115
1116   int nb0d=0, nb2d_tria=0, nb2d_quad=0;
1117   bool isQuadratic = false, ok = true;
1118   if(nbe==1)
1119   {
1120     // C1 must be a circle
1121     ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1122     if(ok) {
1123       const vector<int>& aVec = aResMap[aMesh.GetSubMesh(CircEdge)];
1124       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1125       if(isQuadratic) {
1126         // main nodes
1127         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1128         // radial medium nodes
1129         nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
1130         // other medium nodes
1131         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1132       }
1133       else {
1134         nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1135       }
1136       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1137       nb2d_quad = nb0d;
1138     }
1139   }
1140   else if(nbe==2 && LinEdge1.Orientation() != TopAbs_INTERNAL)
1141   {
1142     // one curve must be a half of circle and other curve must be
1143     // a segment of line
1144     double fp, lp;
1145     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
1146     if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
1147       // not half of circle
1148       return error(COMPERR_BAD_SHAPE);
1149     }
1150     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1151     if( aLine.IsNull() ) {
1152       // other curve not line
1153       return error(COMPERR_BAD_SHAPE);
1154     }
1155     ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1) );
1156     if ( !ok ) {
1157       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1158       ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1159     }
1160     if(ok) {
1161       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1162     }
1163     if(ok) {
1164       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1165       isQuadratic = aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge];
1166       if(isQuadratic) {
1167         // main nodes
1168         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1169         // radial medium nodes
1170         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1171         // other medium nodes
1172         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1173       }
1174       else {
1175         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1176       }
1177       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1178       nb2d_quad = nb2d_tria * myLayerPositions.size();
1179       // add evaluation for edges
1180       vector<int> aResVec(SMDSEntity_Last,0);
1181       if(isQuadratic) {
1182         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
1183         aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
1184       }
1185       else {
1186         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1187         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
1188       }
1189       aResMap[ aMesh.GetSubMesh(LinEdge1) ] = aResVec;
1190     }
1191   }
1192   else  // nbe==3 or ( nbe==2 && linEdge is INTERNAL )
1193   {
1194     if (nbe==2 && LinEdge1.Orientation() == TopAbs_INTERNAL )
1195       LinEdge2 = LinEdge1;
1196
1197     // one curve must be a part of circle and other curves must be
1198     // segments of line
1199     Handle(Geom_Line)  aLine1 = Handle(Geom_Line)::DownCast( getCurve( LinEdge1 ));
1200     Handle(Geom_Line)  aLine2 = Handle(Geom_Line)::DownCast( getCurve( LinEdge2 ));
1201     if( aLine1.IsNull() || aLine2.IsNull() ) {
1202       // other curve not line
1203       return error(COMPERR_BAD_SHAPE);
1204     }
1205     int nbLayers = myLayerPositions.size();
1206     computeLayerPositions( P0, P1, LinEdge2 );
1207     if ( nbLayers != myLayerPositions.size() )
1208       return error("Different hypotheses apply to radial edges");
1209
1210     bool ok = !aResMap.count( aMesh.GetSubMesh(LinEdge1));
1211     if ( !ok ) {
1212       if ( myDistributionHypo || myNbLayerHypo )
1213         ok = true; // override other 1d hyps
1214       else {
1215         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge1) ];
1216         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1217       }
1218     }
1219     if( ok && aResMap.count( aMesh.GetSubMesh(LinEdge2) )) {
1220       if ( myDistributionHypo || myNbLayerHypo )
1221         ok = true; // override other 1d hyps
1222       else {
1223         const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(LinEdge2) ];
1224         ok = ( aVec[SMDSEntity_Node] == myLayerPositions.size() );
1225       }
1226     }
1227     if(ok) {
1228       ok = algo1d->EvaluateCircularEdge( aMesh, CircEdge, aResMap );
1229     }
1230     if(ok) {
1231       const vector<int>& aVec = aResMap[ aMesh.GetSubMesh(CircEdge) ];
1232       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1233       if(isQuadratic) {
1234         // main nodes
1235         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1236         // radial medium nodes
1237         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1238         // other medium nodes
1239         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1240       }
1241       else {
1242         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1243       }
1244       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1245       nb2d_quad = nb2d_tria * myLayerPositions.size();
1246       // add evaluation for edges
1247       vector<int> aResVec(SMDSEntity_Last, 0);
1248       if(isQuadratic) {
1249         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1250         aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1251       }
1252       else {
1253         aResVec[SMDSEntity_Node] = myLayerPositions.size();
1254         aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1255       }
1256       sm = aMesh.GetSubMesh(LinEdge1);
1257       aResMap[sm] = aResVec;
1258       sm = aMesh.GetSubMesh(LinEdge2);
1259       aResMap[sm] = aResVec;
1260     }
1261   }
1262
1263   if(nb0d>0) {
1264     aResVec[0] = nb0d;
1265     if(isQuadratic) {
1266       aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1267       aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1268     }
1269     else {
1270       aResVec[SMDSEntity_Triangle] = nb2d_tria;
1271       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1272     }
1273     return true;
1274   }
1275
1276   // invalid case
1277   sm = aMesh.GetSubMesh(aShape);
1278   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1279   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1280                                         "Submesh can not be evaluated",this));
1281   return false;
1282
1283 }
1284
1285 //================================================================================
1286 /*!
1287  * \brief Return true if applied compute mesh on this shape
1288  */
1289 //================================================================================
1290
1291 bool StdMeshers_RadialQuadrangle_1D2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
1292 {
1293   int nbFoundFaces = 0;
1294   for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces ){
1295     TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
1296     int nbe = analyseFace( exp.Current(), CircEdge, LinEdge1, LinEdge2 );
1297     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge ));
1298     bool ok = ( nbe <= 3 && nbe >= 1 && !aCirc.IsNull() );
1299     if( toCheckAll  && !ok ) return false;
1300     if( !toCheckAll && ok  ) return true;
1301   }
1302   if( toCheckAll && nbFoundFaces != 0 ) return true;
1303   return false;
1304 };