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