src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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