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 <BRepTools.hxx>
  49 #include <BRep_Tool.hxx>
  50 #include <TopExp_Explorer.hxx>
  51 #include <TopoDS.hxx>
  52 //#include <TopoDS_Shell.hxx>
  53 //#include <TopoDS_Solid.hxx>
  54 //#include <TopTools_MapOfShape.hxx>
  55 //#include <gp.hxx>
  56 //#include <gp_Pnt.hxx>
  57
  58
  59 #include <Geom_TrimmedCurve.hxx>
  60 #include <Geom_Circle.hxx>
  61 #include <Geom_Line.hxx>
  62 #include <TColgp_SequenceOfPnt.hxx>
  63 #include <TColgp_SequenceOfPnt2d.hxx>
  64 #include <GeomAPI_ProjectPointOnSurf.hxx>
  65
  66
  67 using namespace std;
  68
  69 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  70 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  71
  72 //typedef StdMeshers_ProjectionUtils TAssocTool;
  73
  74
  75 //=======================================================================
  76 //function : StdMeshers_RadialQuadrangle_1D2D
  77 //purpose  :
  78 //=======================================================================
  79
  80 StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
  81                                                                    int studyId,
  82                                                                    SMESH_Gen* gen)
  83   :SMESH_2D_Algo(hypId, studyId, gen)
  84 {
  85   _name = "RadialQuadrangle_1D2D";
  86   _shapeType = (1 << TopAbs_FACE);        // 1 bit per shape type
  87
  88   _compatibleHypothesis.push_back("LayerDistribution2D");
  89   _compatibleHypothesis.push_back("NumberOfLayers2D");
  90   myNbLayerHypo = 0;
  91   myDistributionHypo = 0;
  92   _requireDescretBoundary = false;
  93 }
  94
  95
  96 //================================================================================
  97 /*!
  98  * \brief Destructor
  99  */
 100 //================================================================================
 101
 102 StdMeshers_RadialQuadrangle_1D2D::~StdMeshers_RadialQuadrangle_1D2D()
 103 {}
 104
 105
 106 //=======================================================================
 107 //function : CheckHypothesis
 108 //purpose  :
 109 //=======================================================================
 110
 111 bool StdMeshers_RadialQuadrangle_1D2D::CheckHypothesis
 112                            (SMESH_Mesh&                          aMesh,
 113                             const TopoDS_Shape&                  aShape,
 114                             SMESH_Hypothesis::Hypothesis_Status& aStatus)
 115 {
 116   // check aShape
 117   myNbLayerHypo = 0;
 118   myDistributionHypo = 0;
 119
 120   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 121
 122   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 123   if ( hyps.size() == 0 ) {
 124     aStatus = SMESH_Hypothesis::HYP_MISSING;
 125     return false;  // can't work with no hypothesis
 126   }
 127
 128   if ( hyps.size() > 1 ) {
 129     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 130     return false;
 131   }
 132
 133   const SMESHDS_Hypothesis *theHyp = hyps.front();
 134
 135   string hypName = theHyp->GetName();
 136
 137   if (hypName == "NumberOfLayers2D") {
 138     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 139     aStatus = SMESH_Hypothesis::HYP_OK;
 140     return true;
 141   }
 142   if (hypName == "LayerDistribution2D") {
 143     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 144     aStatus = SMESH_Hypothesis::HYP_OK;
 145     return true;
 146   }
 147   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 148   return true;
 149 }
 150
 151
 152 //=======================================================================
 153 //function : Compute
 154 //purpose  :
 155 //=======================================================================
 156
 157 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh& aMesh,
 158                                                const TopoDS_Shape& aShape)
 159 {
 160   TopExp_Explorer exp;
 161   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 162
 163   myHelper = new SMESH_MesherHelper( aMesh );
 164   myHelper->IsQuadraticSubMesh( aShape );
 165
 166   myLayerPositions.clear();
 167
 168   TopoDS_Edge E1,E2,E3;
 169   Handle(Geom_Curve) C1,C2,C3;
 170   double f1,l1,f2,l2,f3,l3;
 171   int nbe = 0;
 172   for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 173     nbe++;
 174     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
 175     if(nbe==1) {
 176       E1 = E;
 177       C1 = BRep_Tool::Curve(E,f1,l1);
 178     }
 179     else if(nbe==2) {
 180       E2 = E;
 181       C2 = BRep_Tool::Curve(E,f2,l2);
 182     }
 183     else if(nbe==3) {
 184       E3 = E;
 185       C3 = BRep_Tool::Curve(E,f3,l3);
 186     }
 187   }
 188
 189   if(nbe>3)
 190     return error(COMPERR_BAD_SHAPE);
 191
 192   gp_Pnt P0,P1;
 193   // points for rotation
 194   TColgp_SequenceOfPnt Points;
 195   // angles for rotation
 196   TColStd_SequenceOfReal Angles;
 197   // Nodes1 and Nodes2 - nodes along radiuses
 198   // CNodes - nodes on circle edge
 199   std::vector< const SMDS_MeshNode* > Nodes1, Nodes2, CNodes;
 200   SMDS_MeshNode * NC;
 201   // parameters edge nodes on face
 202   TColgp_SequenceOfPnt2d Pnts2d1, Pnts2d2;
 203   gp_Pnt2d PC;
 204
 205   int faceID = meshDS->ShapeToIndex(aShape);
 206   TopoDS_Face F = TopoDS::Face(aShape);
 207   Handle(Geom_Surface) S = BRep_Tool::Surface(F);
 208
 209   // orientation
 210   bool IsForward = F.Orientation()==TopAbs_FORWARD;
 211
 212   //cout<<"RadialQuadrangle_1D2D::Compute   nbe = "<<nbe<<endl;
 213   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 214   if(nbe==1) {
 215     // C1 must be a circle
 216     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
 217     if( aCirc.IsNull() )
 218       return error(COMPERR_BAD_SHAPE);
 219     CircEdge = E1;
 220
 221     bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
 222     if( !ok ) return false;
 223     std::map< double, const SMDS_MeshNode* > theNodes;
 224     ok = GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 225     if( !ok ) return false;
 226
 227     CNodes.clear();
 228     std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 229     const SMDS_MeshNode* NF = (*itn).second;
 230     CNodes.push_back( (*itn).second );
 231     double fang = (*itn).first;
 232     if ( itn != theNodes.end() ) {
 233       itn++;
 234       for(; itn != theNodes.end(); itn++ ) {
 235         CNodes.push_back( (*itn).second );
 236         double ang = (*itn).first - fang;
 237         if( ang>PI ) ang = ang - 2*PI;
 238         if( ang<-PI ) ang = ang + 2*PI;
 239         Angles.Append( ang );
 240       }
 241     }
 242     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 243     P0 = aCirc->Location();
 244
 245     myLayerPositions.clear();
 246     computeLayerPositions(P0,P1);
 247
 248     exp.Init( CircEdge, TopAbs_VERTEX );
 249     TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
 250     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
 251
 252     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 253     GeomAPI_ProjectPointOnSurf PPS(P0,S);
 254     double U0,V0;
 255     PPS.Parameters(1,U0,V0);
 256     meshDS->SetNodeOnFace(NC, faceID, U0, V0);
 257     PC = gp_Pnt2d(U0,V0);
 258
 259     gp_Vec aVec(P0,P1);
 260     gp_Vec2d aVec2d(PC,p2dV);
 261     Nodes1.resize( myLayerPositions.size()+1 );
 262     Nodes2.resize( myLayerPositions.size()+1 );
 263     int i = 0;
 264     for(; i<myLayerPositions.size(); i++) {
 265       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 266                 P0.Y() + aVec.Y()*myLayerPositions[i],
 267                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 268       Points.Append(P);
 269       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 270       Nodes1[i] = node;
 271       Nodes2[i] = node;
 272       double U = PC.X() + aVec2d.X()*myLayerPositions[i];
 273       double V = PC.Y() + aVec2d.Y()*myLayerPositions[i];
 274       meshDS->SetNodeOnFace( node, faceID, U, V );
 275       Pnts2d1.Append(gp_Pnt2d(U,V));
 276       Pnts2d2.Append(gp_Pnt2d(U,V));
 277     }
 278     Nodes1[Nodes1.size()-1] = NF;
 279     Nodes2[Nodes1.size()-1] = NF;
 280   }
 281   else if(nbe==2) {
 282     // one curve must be a half of circle and other curve must be
 283     // a segment of line
 284     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 285     while( !tc.IsNull() ) {
 286       C1 = tc->BasisCurve();
 287       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 288     }
 289     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 290     while( !tc.IsNull() ) {
 291       C2 = tc->BasisCurve();
 292       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 293     }
 294
 295     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
 296     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
 297     CircEdge = E1;
 298     LinEdge1 = E2;
 299     double fp = f1;
 300     double lp = l1;
 301     if( aCirc.IsNull() ) {
 302       aCirc = Handle(Geom_Circle)::DownCast(C2);
 303       CircEdge = E2;
 304       LinEdge1 = E1;
 305       fp = f2;
 306       lp = l2;
 307       aLine = Handle(Geom_Line)::DownCast(C3);
 308     }
 309     if( aCirc.IsNull() ) {
 310       // not circle
 311       return error(COMPERR_BAD_SHAPE);
 312     }
 313     if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
 314       // not half of circle
 315       return error(COMPERR_BAD_SHAPE);
 316     }
 317     if( aLine.IsNull() ) {
 318       // other curve not line
 319       return error(COMPERR_BAD_SHAPE);
 320     }
 321     SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
 322     if( sm1 ) {
 323       SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS();
 324       if( sdssm1 ) {
 325         if( sm1->GetSubMeshDS()->NbNodes()>0 ) {
 326           SMESH_subMesh* sm = aMesh.GetSubMesh(F);
 327           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 328           smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 329                                                "Invalid set of hypothesises",this));
 330           return false;
 331         }
 332       }
 333     }
 334
 335     bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
 336     if( !ok ) return false;
 337     std::map< double, const SMDS_MeshNode* > theNodes;
 338     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 339
 340     CNodes.clear();
 341     std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 342     const SMDS_MeshNode* NF = (*itn).second;
 343     CNodes.push_back( (*itn).second );
 344     double fang = (*itn).first;
 345     itn++;
 346     const SMDS_MeshNode* NL;
 347     int nbn = 1;
 348     for(; itn != theNodes.end(); itn++ ) {
 349       nbn++;
 350       if( nbn == theNodes.size() )
 351         NL = (*itn).second;
 352       CNodes.push_back( (*itn).second );
 353       double ang = (*itn).first - fang;
 354       if( ang>PI ) ang = ang - 2*PI;
 355       if( ang<-PI ) ang = ang + 2*PI;
 356       Angles.Append( ang );
 357     }
 358     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 359     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 360     P0 = aCirc->Location();
 361
 362     myLayerPositions.clear();
 363     computeLayerPositions(P0,P1);
 364
 365     gp_Vec aVec(P0,P1);
 366     int edgeID = meshDS->ShapeToIndex(LinEdge1);
 367     // check orientation
 368     Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 369     gp_Pnt Ptmp;
 370     Crv->D0(fp,Ptmp);
 371     bool ori = true;
 372     if( P1.Distance(Ptmp) > Precision::Confusion() )
 373       ori = false;
 374     // get UV points for edge
 375     gp_Pnt2d PF,PL;
 376     BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 377     PC = gp_Pnt2d( (PF.X()+PL.X())/2, (PF.Y()+PL.Y())/2 );
 378     gp_Vec2d V2d;
 379     if(ori) V2d = gp_Vec2d(PC,PF);
 380     else V2d = gp_Vec2d(PC,PL);
 381     // add nodes on edge
 382     double cp = (fp+lp)/2;
 383     double dp2 = (lp-fp)/2;
 384     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 385     meshDS->SetNodeOnEdge(NC, edgeID, cp);
 386     Nodes1.resize( myLayerPositions.size()+1 );
 387     Nodes2.resize( myLayerPositions.size()+1 );
 388     int i = 0;
 389     for(; i<myLayerPositions.size(); i++) {
 390       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 391                 P0.Y() + aVec.Y()*myLayerPositions[i],
 392                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 393       Points.Append(P);
 394       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 395       Nodes1[i] = node;
 396       double param;
 397       if(ori)
 398         param = fp + dp2*(1-myLayerPositions[i]);
 399       else
 400         param = cp + dp2*myLayerPositions[i];
 401       meshDS->SetNodeOnEdge(node, edgeID, param);
 402       P = gp_Pnt( P0.X() - aVec.X()*myLayerPositions[i],
 403                   P0.Y() - aVec.Y()*myLayerPositions[i],
 404                   P0.Z() - aVec.Z()*myLayerPositions[i] );
 405       node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 406       Nodes2[i] = node;
 407       if(!ori)
 408         param = fp + dp2*(1-myLayerPositions[i]);
 409       else
 410         param = cp + dp2*myLayerPositions[i];
 411       meshDS->SetNodeOnEdge(node, edgeID, param);
 412       // parameters on face
 413       gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 414                     PC.Y() + V2d.Y()*myLayerPositions[i] );
 415       Pnts2d1.Append(P2d);
 416       P2d = gp_Pnt2d( PC.X() - V2d.X()*myLayerPositions[i],
 417                       PC.Y() - V2d.Y()*myLayerPositions[i] );
 418       Pnts2d2.Append(P2d);
 419     }
 420     Nodes1[ myLayerPositions.size() ] = NF;
 421     Nodes2[ myLayerPositions.size() ] = NL;
 422     // create 1D elements on edge
 423     std::vector< const SMDS_MeshNode* > tmpNodes;
 424     tmpNodes.resize(2*Nodes1.size()+1);
 425     for(i=0; i<Nodes2.size(); i++)
 426       tmpNodes[Nodes2.size()-i-1] = Nodes2[i];
 427     tmpNodes[Nodes2.size()] = NC;
 428     for(i=0; i<Nodes1.size(); i++)
 429       tmpNodes[Nodes2.size()+1+i] = Nodes1[i];
 430     for(i=1; i<tmpNodes.size(); i++) {
 431       SMDS_MeshEdge* ME = myHelper->AddEdge( tmpNodes[i-1], tmpNodes[i] );
 432       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 433     }
 434   }
 435   else { // nbe==3
 436     // one curve must be a part of circle and other curves must be
 437     // segments of line
 438     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 439     while( !tc.IsNull() ) {
 440       C1 = tc->BasisCurve();
 441       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 442     }
 443     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 444     while( !tc.IsNull() ) {
 445       C2 = tc->BasisCurve();
 446       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 447     }
 448     tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
 449     while( !tc.IsNull() ) {
 450       C3 = tc->BasisCurve();
 451       tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
 452     }
 453
 454     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
 455     Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
 456     Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
 457     CircEdge = E1;
 458     LinEdge1 = E2;
 459     LinEdge2 = E3;
 460     double fp = f1;
 461     double lp = l1;
 462     if( aCirc.IsNull() ) {
 463       aCirc = Handle(Geom_Circle)::DownCast(C2);
 464       CircEdge = E2;
 465       LinEdge1 = E3;
 466       LinEdge2 = E1;
 467       fp = f2;
 468       lp = l2;
 469       aLine1 = Handle(Geom_Line)::DownCast(C3);
 470       aLine2 = Handle(Geom_Line)::DownCast(C1);
 471       if( aCirc.IsNull() ) {
 472         aCirc = Handle(Geom_Circle)::DownCast(C3);
 473         CircEdge = E3;
 474         LinEdge1 = E1;
 475         LinEdge2 = E2;
 476         fp = f3;
 477         lp = l3;
 478         aLine1 = Handle(Geom_Line)::DownCast(C1);
 479         aLine2 = Handle(Geom_Line)::DownCast(C2);
 480       }
 481     }
 482     if( aCirc.IsNull() ) {
 483       // not circle
 484       return error(COMPERR_BAD_SHAPE);
 485     }
 486     if( aLine1.IsNull() || aLine2.IsNull() ) {
 487       // other curve not line
 488       return error(COMPERR_BAD_SHAPE);
 489     }
 490     SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
 491     SMESH_subMesh* sm2 = aMesh.GetSubMesh(LinEdge2);
 492     if( sm1 && sm2 ) {
 493       SMESHDS_SubMesh* sdssm1 = sm1->GetSubMeshDS();
 494       SMESHDS_SubMesh* sdssm2 = sm2->GetSubMeshDS();
 495       if( sdssm1 && sdssm2 ) {
 496         if( sm1->GetSubMeshDS()->NbNodes()>0 || sm2->GetSubMeshDS()->NbNodes()>0 ) {
 497           SMESH_subMesh* sm = aMesh.GetSubMesh(F);
 498           SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 499           smError.reset(new SMESH_ComputeError(COMPERR_ALGO_FAILED,
 500                                                "Invalid set of hypothesises",this));
 501           return false;
 502         }
 503       }
 504     }
 505
 506     bool ok = _gen->Compute( aMesh, CircEdge, false, MeshDim_1D );
 507     if( !ok ) return false;
 508     std::map< double, const SMDS_MeshNode* > theNodes;
 509     GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes);
 510
 511     CNodes.clear();
 512     std::map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
 513     const SMDS_MeshNode* NF = (*itn).second;
 514     CNodes.push_back( (*itn).second );
 515     double fang = (*itn).first;
 516     itn++;
 517     const SMDS_MeshNode* NL;
 518     int nbn = 1;
 519     for(; itn != theNodes.end(); itn++ ) {
 520       nbn++;
 521       if( nbn == theNodes.size() )
 522         NL = (*itn).second;
 523       CNodes.push_back( (*itn).second );
 524       double ang = (*itn).first - fang;
 525       if( ang>PI ) ang = ang - 2*PI;
 526       if( ang<-PI ) ang = ang + 2*PI;
 527       Angles.Append( ang );
 528     }
 529     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
 530     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
 531     P0 = aCirc->Location();
 532
 533     myLayerPositions.clear();
 534     computeLayerPositions(P0,P1);
 535
 536     exp.Init( LinEdge1, TopAbs_VERTEX );
 537     TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
 538     exp.Next();
 539     TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
 540     gp_Pnt PE1 = BRep_Tool::Pnt(V1);
 541     gp_Pnt PE2 = BRep_Tool::Pnt(V2);
 542     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 543         ( P1.Distance(PE2) > Precision::Confusion() ) ) {
 544       TopoDS_Edge E = LinEdge1;
 545       LinEdge1 = LinEdge2;
 546       LinEdge2 = E;
 547     }
 548     TopoDS_Vertex VC;
 549     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
 550         ( P2.Distance(PE1) > Precision::Confusion() ) ) {
 551       VC = V1;
 552     }
 553     else VC = V2;
 554     int vertID = meshDS->ShapeToIndex(VC);
 555     // LinEdge1
 556     int edgeID = meshDS->ShapeToIndex(LinEdge1);
 557     gp_Vec aVec(P0,P1);
 558     // check orientation
 559     Handle(Geom_Curve) Crv = BRep_Tool::Curve(LinEdge1,fp,lp);
 560     gp_Pnt Ptmp;
 561     Crv->D0(fp,Ptmp);
 562     bool ori = false;
 563     if( P1.Distance(Ptmp) > Precision::Confusion() )
 564       ori = true;
 565     // get UV points for edge
 566     gp_Pnt2d PF,PL;
 567     BRep_Tool::UVPoints( LinEdge1, TopoDS::Face(aShape), PF, PL );
 568     gp_Vec2d V2d;
 569     if(ori) {
 570       V2d = gp_Vec2d(PF,PL);
 571       PC = PF;
 572     }
 573     else {
 574       V2d = gp_Vec2d(PL,PF);
 575       PC = PL;
 576     }
 577     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
 578     meshDS->SetNodeOnVertex(NC, vertID);
 579     double dp = lp-fp;
 580     Nodes1.resize( myLayerPositions.size()+1 );
 581     int i = 0;
 582     for(; i<myLayerPositions.size(); i++) {
 583       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 584                 P0.Y() + aVec.Y()*myLayerPositions[i],
 585                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 586       Points.Append(P);
 587       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 588       Nodes1[i] = node;
 589       double param;
 590       if(!ori)
 591         param = fp + dp*(1-myLayerPositions[i]);
 592       else
 593         param = fp + dp*myLayerPositions[i];
 594       meshDS->SetNodeOnEdge(node, edgeID, param);
 595       // parameters on face
 596       gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 597                     PC.Y() + V2d.Y()*myLayerPositions[i] );
 598       Pnts2d1.Append(P2d);
 599     }
 600     Nodes1[ myLayerPositions.size() ] = NF;
 601     // create 1D elements on edge
 602     SMDS_MeshEdge* ME = myHelper->AddEdge( NC, Nodes1[0] );
 603     if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 604     for(i=1; i<Nodes1.size(); i++) {
 605       SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes1[i-1], Nodes1[i] );
 606       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 607     }
 608     // LinEdge2
 609     edgeID = meshDS->ShapeToIndex(LinEdge2);
 610     aVec = gp_Vec(P0,P2);
 611     // check orientation
 612     Crv = BRep_Tool::Curve(LinEdge2,fp,lp);
 613     Crv->D0(fp,Ptmp);
 614     ori = false;
 615     if( P2.Distance(Ptmp) > Precision::Confusion() )
 616       ori = true;
 617     // get UV points for edge
 618     BRep_Tool::UVPoints( LinEdge2, TopoDS::Face(aShape), PF, PL );
 619     if(ori) {
 620       V2d = gp_Vec2d(PF,PL);
 621       PC = PF;
 622     }
 623     else {
 624       V2d = gp_Vec2d(PL,PF);
 625       PC = PL;
 626     }
 627     dp = lp-fp;
 628     Nodes2.resize( myLayerPositions.size()+1 );
 629     for(i=0; i<myLayerPositions.size(); i++) {
 630       gp_Pnt P( P0.X() + aVec.X()*myLayerPositions[i],
 631                 P0.Y() + aVec.Y()*myLayerPositions[i],
 632                 P0.Z() + aVec.Z()*myLayerPositions[i] );
 633       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
 634       Nodes2[i] = node;
 635       double param;
 636       if(!ori)
 637         param = fp + dp*(1-myLayerPositions[i]);
 638       else
 639         param = fp + dp*myLayerPositions[i];
 640       meshDS->SetNodeOnEdge(node, edgeID, param);
 641       // parameters on face
 642       gp_Pnt2d P2d( PC.X() + V2d.X()*myLayerPositions[i],
 643                     PC.Y() + V2d.Y()*myLayerPositions[i] );
 644       Pnts2d2.Append(P2d);
 645     }
 646     Nodes2[ myLayerPositions.size() ] = NL;
 647     // create 1D elements on edge
 648     ME = myHelper->AddEdge( NC, Nodes2[0] );
 649     if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 650     for(i=1; i<Nodes2.size(); i++) {
 651       SMDS_MeshEdge* ME = myHelper->AddEdge( Nodes2[i-1], Nodes2[i] );
 652       if(ME) meshDS->SetMeshElementOnShape(ME, edgeID);
 653     }
 654   }
 655
 656   // create nodes and mesh elements on face
 657   // find axis of rotation
 658   gp_Pnt P2 = gp_Pnt( CNodes[1]->X(), CNodes[1]->Y(), CNodes[1]->Z() );
 659   gp_Vec Vec1(P0,P1);
 660   gp_Vec Vec2(P0,P2);
 661   gp_Vec Axis = Vec1.Crossed(Vec2);
 662   // create elements
 663   int i = 1;
 664   //cout<<"Angles.Length() = "<<Angles.Length()<<"   Points.Length() = "<<Points.Length()<<endl;
 665   //cout<<"Nodes1.size() = "<<Nodes1.size()<<"   Pnts2d1.Length() = "<<Pnts2d1.Length()<<endl;
 666   for(; i<Angles.Length(); i++) {
 667     std::vector< const SMDS_MeshNode* > tmpNodes;
 668     tmpNodes.reserve(Nodes1.size());
 669     gp_Trsf aTrsf;
 670     gp_Ax1 theAxis(P0,gp_Dir(Axis));
 671     aTrsf.SetRotation( theAxis, Angles.Value(i) );
 672     gp_Trsf2d aTrsf2d;
 673     aTrsf2d.SetRotation( PC, Angles.Value(i) );
 674     // create nodes
 675     int j = 1;
 676     for(; j<=Points.Length(); j++) {
 677       double cx,cy,cz;
 678       Points.Value(j).Coord( cx, cy, cz );
 679       aTrsf.Transforms( cx, cy, cz );
 680       SMDS_MeshNode* node = myHelper->AddNode( cx, cy, cz );
 681       // find parameters on face
 682       Pnts2d1.Value(j).Coord( cx, cy );
 683       aTrsf2d.Transforms( cx, cy );
 684       // set node on face
 685       meshDS->SetNodeOnFace( node, faceID, cx, cy );
 686       tmpNodes[j-1] = node;
 687     }
 688     // create faces
 689     tmpNodes[Points.Length()] = CNodes[i];
 690     // quad
 691     for(j=0; j<Nodes1.size()-1; j++) {
 692       SMDS_MeshFace* MF;
 693       if(IsForward)
 694         MF = myHelper->AddFace( tmpNodes[j], Nodes1[j],
 695                                 Nodes1[j+1], tmpNodes[j+1] );
 696       else
 697         MF = myHelper->AddFace( tmpNodes[j], tmpNodes[j+1],
 698                                 Nodes1[j+1], Nodes1[j] );
 699       if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 700     }
 701     // tria
 702     SMDS_MeshFace* MF;
 703     if(IsForward)
 704       MF = myHelper->AddFace( NC, Nodes1[0], tmpNodes[0] );
 705     else
 706       MF = myHelper->AddFace( NC, tmpNodes[0], Nodes1[0] );
 707     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 708     for(j=0; j<Nodes1.size(); j++) {
 709       Nodes1[j] = tmpNodes[j];
 710     }
 711   }
 712   // create last faces
 713   // quad
 714   for(i=0; i<Nodes1.size()-1; i++) {
 715     SMDS_MeshFace* MF;
 716     if(IsForward)
 717       MF = myHelper->AddFace( Nodes2[i], Nodes1[i],
 718                               Nodes1[i+1], Nodes2[i+1] );
 719     else
 720       MF = myHelper->AddFace( Nodes2[i],  Nodes2[i+1],
 721                               Nodes1[i+1], Nodes1[i] );
 722     if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 723   }
 724   // tria
 725   SMDS_MeshFace* MF;
 726   if(IsForward)
 727     MF = myHelper->AddFace( NC, Nodes1[0], Nodes2[0] );
 728   else
 729     MF = myHelper->AddFace( NC, Nodes2[0], Nodes1[0] );
 730   if(MF) meshDS->SetMeshElementOnShape(MF, faceID);
 731
 732
 733   // to delete helper at exit from Compute()
 734   std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
 735
 736   return true;
 737 }
 738
 739
 740 //================================================================================
 741 //================================================================================
 742 /*!
 743  * \brief Class computing layers distribution using data of
 744  *        StdMeshers_LayerDistribution hypothesis
 745  */
 746 //================================================================================
 747 //================================================================================
 748
 749 class TNodeDistributor: public StdMeshers_Regular_1D
 750 {
 751   list <const SMESHDS_Hypothesis *> myUsedHyps;
 752 public:
 753   // -----------------------------------------------------------------------------
 754   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 755   {
 756     const int myID = -1000;
 757     map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
 758     map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
 759     if ( id_algo == algoMap.end() )
 760       return new TNodeDistributor( myID, 0, aMesh.GetGen() );
 761     return static_cast< TNodeDistributor* >( id_algo->second );
 762   }
 763   // -----------------------------------------------------------------------------
 764   bool Compute( vector< double > &                  positions,
 765                 gp_Pnt                              pIn,
 766                 gp_Pnt                              pOut,
 767                 SMESH_Mesh&                         aMesh,
 768                 const StdMeshers_LayerDistribution* hyp)
 769   {
 770     double len = pIn.Distance( pOut );
 771     if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
 772
 773     if ( !hyp || !hyp->GetLayerDistribution() )
 774       return error( "Invalid LayerDistribution hypothesis");
 775     myUsedHyps.clear();
 776     myUsedHyps.push_back( hyp->GetLayerDistribution() );
 777
 778     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
 779     SMESH_Hypothesis::Hypothesis_Status aStatus;
 780     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 781       return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
 782                     "with LayerDistribution hypothesis");
 783
 784     BRepAdaptor_Curve C3D(edge);
 785     double f = C3D.FirstParameter(), l = C3D.LastParameter();
 786     list< double > params;
 787     if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
 788       return error("StdMeshers_Regular_1D failed to compute layers distribution");
 789
 790     positions.clear();
 791     positions.reserve( params.size() );
 792     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 793       positions.push_back( *itU / len );
 794     return true;
 795   }
 796 protected:
 797   // -----------------------------------------------------------------------------
 798   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
 799     : StdMeshers_Regular_1D( hypId, studyId, gen)
 800   {
 801   }
 802   // -----------------------------------------------------------------------------
 803   virtual const list <const SMESHDS_Hypothesis *> &
 804     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 805   {
 806     return myUsedHyps;
 807   }
 808   // -----------------------------------------------------------------------------
 809 };
 810
 811 //================================================================================
 812 /*!
 813  * \brief Compute positions of nodes between the internal and the external surfaces
 814   * \retval bool - is a success
 815  */
 816 //================================================================================
 817
 818 bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt& pIn,
 819                                                              const gp_Pnt& pOut)
 820 {
 821   if ( myNbLayerHypo )
 822   {
 823     int nbSegments = myNbLayerHypo->GetNumberOfLayers();
 824     myLayerPositions.resize( nbSegments - 1 );
 825     for ( int z = 1; z < nbSegments; ++z )
 826       myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
 827     return true;
 828   }
 829   if ( myDistributionHypo ) {
 830     SMESH_Mesh * mesh = myHelper->GetMesh();
 831     if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
 832                                                             *mesh, myDistributionHypo ))
 833     {
 834       error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
 835       return false;
 836     }
 837   }
 838   RETURN_BAD_RESULT("Bad hypothesis");
 839 }
 840
 841
 842 //=======================================================================
 843 //function : Evaluate
 844 //purpose  :
 845 //=======================================================================
 846
 847 bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
 848                                                 const TopoDS_Shape& aShape,
 849                                                 MapShapeNbElems& aResMap)
 850 {
 851   if( aShape.ShapeType() != TopAbs_FACE ) {
 852     return false;
 853   }
 854   SMESH_subMesh * smf = aMesh.GetSubMesh(aShape);
 855   MapShapeNbElemsItr anIt = aResMap.find(smf);
 856   if( anIt != aResMap.end() ) {
 857     return false;
 858   }
 859
 860   myLayerPositions.clear();
 861   gp_Pnt P0(0,0,0);
 862   gp_Pnt P1(100,0,0);
 863   computeLayerPositions(P0,P1);
 864
 865   TopoDS_Edge E1,E2,E3;
 866   Handle(Geom_Curve) C1,C2,C3;
 867   double f1,l1,f2,l2,f3,l3;
 868   int nbe = 0;
 869   TopExp_Explorer exp;
 870   for ( exp.Init( aShape, TopAbs_EDGE ); exp.More(); exp.Next() ) {
 871     nbe++;
 872     TopoDS_Edge E = TopoDS::Edge( exp.Current() );
 873     if(nbe==1) {
 874       E1 = E;
 875       C1 = BRep_Tool::Curve(E,f1,l1);
 876     }
 877     else if(nbe==2) {
 878       E2 = E;
 879       C2 = BRep_Tool::Curve(E,f2,l2);
 880     }
 881     else if(nbe==3) {
 882       E3 = E;
 883       C3 = BRep_Tool::Curve(E,f3,l3);
 884     }
 885   }
 886
 887   TopoDS_Edge CircEdge, LinEdge1, LinEdge2;
 888   int nb0d=0, nb2d_tria=0, nb2d_quad=0;
 889   bool isQuadratic = false;
 890   if(nbe==1) {
 891     // C1 must be a circle
 892     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
 893     if( !aCirc.IsNull() ) {
 894       bool ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
 895       if(ok) {
 896         SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
 897         MapShapeNbElemsItr anIt = aResMap.find(sm);
 898         std::vector<int> aVec = (*anIt).second;
 899         isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
 900         if(isQuadratic) {
 901           // main nodes
 902           nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
 903           // radial medium nodes
 904           nb0d += (aVec[SMDSEntity_Node]+1) * (myLayerPositions.size()+1);
 905           // other medium nodes
 906           nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
 907         }
 908         else {
 909           nb0d = (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
 910         }
 911         nb2d_tria = aVec[SMDSEntity_Node] + 1;
 912         nb2d_quad = nb0d;
 913       }
 914     }
 915   }
 916   else if(nbe==2) {
 917     // one curve must be a half of circle and other curve must be
 918     // a segment of line
 919     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 920     while( !tc.IsNull() ) {
 921       C1 = tc->BasisCurve();
 922       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 923     }
 924     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 925     while( !tc.IsNull() ) {
 926       C2 = tc->BasisCurve();
 927       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 928     }
 929     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
 930     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(C2);
 931     CircEdge = E1;
 932     LinEdge1 = E2;
 933     double fp = f1;
 934     double lp = l1;
 935     if( aCirc.IsNull() ) {
 936       aCirc = Handle(Geom_Circle)::DownCast(C2);
 937       CircEdge = E2;
 938       LinEdge1 = E1;
 939       fp = f2;
 940       lp = l2;
 941       aLine = Handle(Geom_Line)::DownCast(C3);
 942     }
 943     bool ok = !aCirc.IsNull() && !aLine.IsNull();
 944     if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
 945       // not half of circle
 946       ok = false;
 947     }
 948     SMESH_subMesh* sm1 = aMesh.GetSubMesh(LinEdge1);
 949     MapShapeNbElemsItr anIt = aResMap.find(sm1);
 950     if( anIt!=aResMap.end() ) {
 951       ok = false;
 952     }
 953     if(ok) {
 954       ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
 955     }
 956     if(ok) {
 957       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
 958       MapShapeNbElemsItr anIt = aResMap.find(sm);
 959       std::vector<int> aVec = (*anIt).second;
 960       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
 961       if(isQuadratic) {
 962         // main nodes
 963         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
 964         // radial medium nodes
 965         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
 966         // other medium nodes
 967         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
 968       }
 969       else {
 970         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
 971       }
 972       nb2d_tria = aVec[SMDSEntity_Node] + 1;
 973       nb2d_quad = nb2d_tria * myLayerPositions.size();
 974       // add evaluation for edges
 975       std::vector<int> aResVec(SMDSEntity_Last);
 976       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 977       if(isQuadratic) {
 978         aResVec[SMDSEntity_Node] = 4*myLayerPositions.size() + 3;
 979         aResVec[SMDSEntity_Quad_Edge] = 2*myLayerPositions.size() + 2;
 980       }
 981       else {
 982         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
 983         aResVec[SMDSEntity_Edge] = 2*myLayerPositions.size() + 2;
 984       }
 985       sm = aMesh.GetSubMesh(LinEdge1);
 986       aResMap.insert(std::make_pair(sm,aResVec));
 987     }
 988   }
 989   else { // nbe==3
 990     // one curve must be a part of circle and other curves must be
 991     // segments of line
 992     Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 993     while( !tc.IsNull() ) {
 994       C1 = tc->BasisCurve();
 995       tc = Handle(Geom_TrimmedCurve)::DownCast(C1);
 996     }
 997     tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
 998     while( !tc.IsNull() ) {
 999       C2 = tc->BasisCurve();
1000       tc = Handle(Geom_TrimmedCurve)::DownCast(C2);
1001     }
1002     tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1003     while( !tc.IsNull() ) {
1004       C3 = tc->BasisCurve();
1005       tc = Handle(Geom_TrimmedCurve)::DownCast(C3);
1006     }
1007     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast(C1);
1008     Handle(Geom_Line) aLine1 = Handle(Geom_Line)::DownCast(C2);
1009     Handle(Geom_Line) aLine2 = Handle(Geom_Line)::DownCast(C3);
1010     CircEdge = E1;
1011     LinEdge1 = E2;
1012     LinEdge2 = E3;
1013     double fp = f1;
1014     double lp = l1;
1015     if( aCirc.IsNull() ) {
1016       aCirc = Handle(Geom_Circle)::DownCast(C2);
1017       CircEdge = E2;
1018       LinEdge1 = E3;
1019       LinEdge2 = E1;
1020       fp = f2;
1021       lp = l2;
1022       aLine1 = Handle(Geom_Line)::DownCast(C3);
1023       aLine2 = Handle(Geom_Line)::DownCast(C1);
1024       if( aCirc.IsNull() ) {
1025         aCirc = Handle(Geom_Circle)::DownCast(C3);
1026         CircEdge = E3;
1027         LinEdge1 = E1;
1028         LinEdge2 = E2;
1029         fp = f3;
1030         lp = l3;
1031         aLine1 = Handle(Geom_Line)::DownCast(C1);
1032         aLine2 = Handle(Geom_Line)::DownCast(C2);
1033       }
1034     }
1035     bool ok = !aCirc.IsNull() && !aLine1.IsNull() && !aLine1.IsNull();
1036     SMESH_subMesh* sm = aMesh.GetSubMesh(LinEdge1);
1037     MapShapeNbElemsItr anIt = aResMap.find(sm);
1038     if( anIt!=aResMap.end() ) {
1039       ok = false;
1040     }
1041     sm = aMesh.GetSubMesh(LinEdge2);
1042     anIt = aResMap.find(sm);
1043     if( anIt!=aResMap.end() ) {
1044       ok = false;
1045     }
1046     if(ok) {
1047       ok = _gen->Evaluate( aMesh, CircEdge, aResMap );
1048     }
1049     if(ok) {
1050       SMESH_subMesh * sm = aMesh.GetSubMesh(CircEdge);
1051       MapShapeNbElemsItr anIt = aResMap.find(sm);
1052       std::vector<int> aVec = (*anIt).second;
1053       isQuadratic = aVec[SMDSEntity_Quad_Edge]>aVec[SMDSEntity_Edge];
1054       if(isQuadratic) {
1055         // main nodes
1056         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1057         // radial medium nodes
1058         nb0d += aVec[SMDSEntity_Node] * (myLayerPositions.size()+1);
1059         // other medium nodes
1060         nb0d += (aVec[SMDSEntity_Node]+1) * myLayerPositions.size();
1061       }
1062       else {
1063         nb0d = aVec[SMDSEntity_Node] * myLayerPositions.size();
1064       }
1065       nb2d_tria = aVec[SMDSEntity_Node] + 1;
1066       nb2d_quad = nb2d_tria * myLayerPositions.size();
1067       // add evaluation for edges
1068       std::vector<int> aResVec(SMDSEntity_Last);
1069       for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1070       if(isQuadratic) {
1071         aResVec[SMDSEntity_Node] = 2*myLayerPositions.size() + 1;
1072         aResVec[SMDSEntity_Quad_Edge] = myLayerPositions.size() + 1;
1073       }
1074       else {
1075         aResVec[SMDSEntity_Node] = myLayerPositions.size();
1076         aResVec[SMDSEntity_Edge] = myLayerPositions.size() + 1;
1077       }
1078       sm = aMesh.GetSubMesh(LinEdge1);
1079       aResMap.insert(std::make_pair(sm,aResVec));
1080       sm = aMesh.GetSubMesh(LinEdge2);
1081       aResMap.insert(std::make_pair(sm,aResVec));
1082     }
1083   }
1084
1085   std::vector<int> aResVec(SMDSEntity_Last);
1086   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1087   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
1088
1089   //cout<<"nb0d = "<<nb0d<<"   nb2d_tria = "<<nb2d_tria<<"   nb2d_quad = "<<nb2d_quad<<endl;
1090   if(nb0d>0) {
1091     aResVec[0] = nb0d;
1092     if(isQuadratic) {
1093       aResVec[SMDSEntity_Quad_Triangle] = nb2d_tria;
1094       aResVec[SMDSEntity_Quad_Quadrangle] = nb2d_quad;
1095     }
1096     else {
1097       aResVec[SMDSEntity_Triangle] = nb2d_tria;
1098       aResVec[SMDSEntity_Quadrangle] = nb2d_quad;
1099     }
1100     aResMap.insert(std::make_pair(sm,aResVec));
1101     return true;
1102   }
1103
1104   // invalid case
1105   aResMap.insert(std::make_pair(sm,aResVec));
1106   SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1107   smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,
1108                                         "Submesh can not be evaluated",this));
1109   return false;
1110
1111 }
1112