src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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