src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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