src/StdMeshers/StdMeshers_RadialQuadrangle_1D2D.cxx

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