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