src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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